# Axion-philic cosmological moduli

Kwang Sik Jeong<EMAIL_ADDRESS>Department of Physics, Pusan National
University, Busan 46241, Korea  Wan Il Park<EMAIL_ADDRESS>Division of
Science Education and Institute of Fusion Science, Jeon-buk National
University, Jeonju, Jeonbuk 54896, Korea

###### Abstract

We show that string moduli have axion-philic nature owing to the model-
insensitive derivative interactions arising from the Kähler potential. The
decay of a modulus into stringy axions occurs without suppression by the mass
of final states. Interestingly, it turns out to hold in general not only for
the scalar partner of the stringy axion but also for any other moduli. The
decay into (pseudo) Nambu-Goldstone bosons (NGBs) also avoids such mass
suppression if the modulus is lighter than or similar in mass to the scalar
partner of the NGB. Such axion-philic nature makes string moduli a natural
source of an observable amount of dark radiation in string compactifications
involving ultralight stringy axions, and possibly in extensions of the
Standard Model that include a cosmologically stable NGB such as the QCD axion.
In the latter case, the fermionic superpartner of the NGB can also contribute
to the dark matter as a feebly interacting massive particle.

††preprint: PNUTP-21-A14

##  I Introduction

String compactification involves many moduli with various masses. Their
Planck-suppressed couplings make them long-lived, and typically the lightest
one most relevant to cosmology. Displaced from the post-inflationary potential
minimum during inflation, generically, moduli are coherently generated after
the end of inflation and come to dominate the energy density of the universe
before their decay. The standard hot thermal universe for a successful Big
Bang Nucleosynthesis (BBN) then requires the lightest modulus to decay
dominantly into the standard model (SM) particles well before the BBN epoch.
Accordingly, the modulus has to be heavier than about $100$ TeV
Coughlan:1983ci ; Ellis:1986zt ; Goncharov:1984qm , and its non-SM decay
channels should be suppressed at least by a couple of orders of magnitude
relative to the SM ones.

As the minimal interaction of a modulus with matter fields, the coupling to
matter kinetic terms is incapable of rendering the modulus to dominantly decay
into the SM sector, questioning the compatibility of string moduli with a
successful cosmology. This is because the modulus decay into matter particles
is generally suppressed by the mass of final states, whereas the decay into
gravitinos or stringy axions has no such suppression. The former may cause the
moduli-induced gravitino problem Endo:2006zj ; Nakamura:2006uc ; Dine:2006ii ,
while the latter may lead to too much dark radiation Cicoli:2012aq ;
Higaki:2012ar ; Higaki:2013lra . A viable reheating after modulus domination
thus necessitates specific interactions between the modulus and the SM sector.
In the Kachru-Kallosh-Linde-Trivedi (KKLT) Kachru:2003aw  and LARGE volume
scenarios Balasubramanian:2005zx , which are the explicit realizations of
four-dimensional de Sitter vacua with all string moduli stabilized, one may
rely on modulus-dependent gauge couplings Nakamura:2008ey  or a modulus
coupling to the Higgs bilinear operator in the Kähler potential Cicoli:2012aq
; Higaki:2012ar ; Higaki:2013lra .

In this paper we point out that the decay of a modulus into axions, which is
induced by the derivative interactions arising from its coupling to the Kähler
potential, can occur without suppression by the mass of final states. This is
the case for the decay into stringy axions regardless of from which modulus
the stringy axion comes. The decay into Nambu-Goldstone bosons (NGBs) also
avoids mass suppression if the modulus is not much heavier than the scalar
partner of the NGB, i.e. the saxion. All these features are insensitive to the
details of moduli stabilization and string compactification, indicating that
string moduli are basically axion-philic. Such axion-philic nature makes a
crucial impact on cosmology because a sizable amount of dark radiation is
produced from their decay in string compactifications giving ultralight
stringy axions, and possibly in SM extensions including a cosmologically
stable axion, like the QCD axion solving the strong CP problem Peccei:1977hh
or a light hidden photon.

Recovering the standard hot big bang universe for the successful BBN requires
additional couplings to the SM sector if the universe undergoes a modulus-
dominated phase. As mentioned, modulus-dependent gauge couplings and/or a
modulus coupling to the holomorphic bilinear operator of SM-charged matter
superfields in the Kähler potential may be necessary. On the other hand, if
the modulus couples to a holomorphic bilinear operator involving the NGB
superfield, NGBs are produced from the modulus decay at the rate without mass
suppression. Also, the fermionic partner of the NGB can contribute to dark
matter as a feebly interacting massive particle (FIMP). If sufficiently light,
it may resolve the moduli-induced gravitino problem.

This paper is organized as follows. In section II, we explore model
independent and dependent interactions of a string modulus inducing its decay
without mass suppression, and show that moduli are generally axion-philic
insensitive to the details of model. The cosmological implications of the
axion-philic nature of moduli are discussed in section III. Section IV is the
conclusions.

##  II Moduli interactions

String moduli generally couple to a (pseudo) NGB and stringy axion, if exists,
through derivative interactions irrespective of the details of string
compactification. We examine how significantly such interactions contribute to
their decay. We also explore model-dependent derivative interactions, and
discuss how to make moduli mainly decay into the SM sector as required for
successful cosmology.

###  II.1 Model-independent axionic decay modes

It is well known that, for a string modulus stabilized while respecting the
associated shift symmetry, the decay of its real component into a pair of its
imaginary one, dubbed as stringy axion, is not mass suppressed Cicoli:2012aq ;
Higaki:2012ar ; Higaki:2013lra . It is also known that the decays of moduli
into matter particles caused by derivative interactions arising from the
Kähler potential are however typically suppressed by the mass of final states.
In this work we show that, interestingly, such mass suppression does not occur
for the decay into NGBs too if the scalar partner of the NGB, i.e. the saxion,
is heavier than or similar in mass to the modulus. We further find that,
regardless of from which modulus the stringy axion originates, moduli decay
into stringy axions without mass suppression.

In order to see the effect of derivative interactions, we consider a simple
model with a single modulus $S$ that enjoys a shift symmetry, ${\rm
Im}(S)\to{\rm Im}(S)+{\rm constant}$, at the perturbative level. The shift
symmetry is not essential for determining the modulus decay width, but makes
our discussion simpler because it ensures that the modulus appears in the
Kähler potential in the combination of $S+S^{*}$. The Kähler potential is
written

$K=K_{0}+Z\Phi\Phi^{*},$ (1)

where $\Phi=\phi+\sqrt{2}\theta\psi+\theta\theta F$ represents a matter
superfield, and we set the reduced Planck mass, $M_{Pl}=1$, unless explicitly
written. Here $K_{0}$ and $Z$ are a function of $S+S^{*}$. Expanding the
modulus around its vacuum expectation value (VEV)

$\delta S\equiv S-S_{0}=s+\sqrt{2}\theta\tilde{s}+\theta\theta F^{S},$ (2)

one always finds the Kähler potential interaction of the modulus with matter
superfields

$(\delta S+\delta S^{*})\Phi\Phi^{*},$ (3)

regardless of the details of the model, where the overall coupling constant is
proportional to $\langle\partial_{S}Z\rangle$. If matter fields do not develop
a VEV, the modulus decay is determined simply by the interactions arising from
the $D$-term of (3). In terms of the component fields, they read

$\left.\delta S\Phi\Phi^{*}\right|_{D}=sFF^{*}+\phi
F^{S}F^{*}-s\phi\partial^{2}\phi^{*}-is\psi\sigma^{\mu}\partial_{\mu}\bar{\psi},$
(4)

omitting total derivatives. Written in the above form, the derivative
interactions are directly combined with the equations of motion for $\phi$ and
$\psi$ to show that the decay rate of the modulus is suppressed by the mass of
decay products.

On the other hand, if $\phi$ develops a nonzero VEV, $\phi$ and $\psi$ are not
the mass eigenstates anymore, and there also arises kinetic mixing between the
moduli and matter fields. To see how the modulus decay is affected, we
consider the case where $\phi$ spontaneously breaks a global U$(1)$ symmetry.
Then the matter scalar field is decomposed into the saxion and the NGB

$\phi(x)=\frac{f}{\sqrt{2}}\left(1+\frac{\sigma(x)}{f}\right)e^{i\frac{a(x)}{f}},$
(5)

with $f$ being the axion decay constant, while the modulus field is written

$s(x)=\frac{1}{\sqrt{2}}\left(s_{r}(x)+is_{i}(x)\right),$ (6)

for the real scalars $s_{r}$ and $s_{i}$. The derivative interaction from
$\delta S\Phi\Phi^{*}$,

$\Delta{\cal L}=\sqrt{2}\kappa\left(s\phi\partial^{2}\phi^{*}+{\rm
h.c.}\right),$ (7)

gives rise to the kinetic mixing terms and scalar cubic interactions

$\displaystyle\kappa^{-1}\Delta{\cal L}|_{2}$ $\displaystyle=$ $\displaystyle
fs_{r}\partial^{2}\sigma+fs_{i}\partial^{2}a,$ (8)
$\displaystyle\kappa^{-1}\Delta{\cal L}|_{3}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}aa\partial^{2}s_{r}+a\sigma\partial^{2}s_{i}$ (9)
$\displaystyle+\,(s_{r}\sigma-
s_{i}a)\partial^{2}\sigma+(s_{r}a+s_{i}\sigma)\partial^{2}a,$

up to total derivatives. Here the overall coupling constant is determined by

$\kappa=-\frac{1}{\sqrt{2}}\left\langle(\partial^{2}_{S}K_{0})^{-\frac{1}{2}}\partial_{S}\ln
Z\right\rangle,$ (10)

and it is generally of the order unity. 111 For Kähler moduli, $\kappa$ is
fixed by the location of the corresponding matter field in extra dimensions.
It is then possible to suppress $\kappa$, but up to by a loop factor because
it receives quantum corrections from non-perturbative effects, string loops,
and higher order $\alpha^{\prime}$ corrections.  As will be shown shortly, the
kinetic mixing has critical impacts on the modulus decay when combined with
the cubic interactions arising from the matter Kähler potential,

$\Phi\Phi^{*}|_{D,3}=\frac{1}{2f}aa\partial^{2}\sigma-\frac{1}{f}\sigma
a\partial^{2}a,$ (11)

up to total derivatives.

The mixing terms (8) should be removed to examine how much the derivative
interactions (9) contribute to the modulus decay. The kinetic mixing between
the CP-odd scalar bosons is eliminated by taking the field transformation

$s_{i}\to s_{i},\quad a\to a+\kappa fs_{i}.$ (12)

After the field transformation, i.e. in the canonical mass basis, all the
derivative cubic couplings of $s_{i}$ vanish because the contributions from
(9) and (11) exactly cancel each other. Hence there remain no interactions
inducing the decay of $s_{i}$. For the CP-even scalar bosons, one can remove
their kinetic mixing by taking the field transformation

$s_{r}\to s_{r}+\epsilon_{1}\kappa
f\sigma,\quad\sigma\to\sigma+\epsilon_{2}\kappa fs_{r},$ (13)

for $f\ll M_{Pl}$, with the coefficients determined by

$\epsilon_{1}\simeq-\frac{m^{2}_{\sigma}}{m^{2}_{s_{r}}-m^{2}_{\sigma}},\quad\epsilon_{2}\simeq\frac{m^{2}_{s_{r}}}{m^{2}_{s_{r}}-m^{2}_{\sigma}},$
(14)

where $m_{i}$ denotes the mass of the indicated scalar boson, and we have
neglected small corrections suppressed by $(f/M_{Pl})^{2}$. As a result, from
(9) and (11), the derivative interactions for $s_{r}$ read

$\Delta{\cal L}|_{3}=\kappa
s_{r}\sigma\partial^{2}\sigma+\frac{\kappa}{2}\frac{m^{2}_{\sigma}}{m^{2}_{s_{r}}-m^{2}_{\sigma}}aa\partial^{2}s_{r},$
(15)

in the canonical mass basis. Therefore, if $m_{s_{r}}\gg m_{\sigma}$, the
equations of motion lead to that both $s_{r}\to aa$ and $s_{r}\to\sigma\sigma$
are suppressed by the saxion mass. On the other hand, if $m_{s_{r}}\ll
m_{\sigma}$, i.e. if the saxion is heavier than the modulus, the relevant
interaction for the modulus decay becomes 222 Eq. (14) is valid for
$f/M_{Pl}\ll|m^{2}_{s_{r}}-m^{2}_{\sigma}|/(m^{2}_{s_{r}}+m^{2}_{\sigma})$.
Although it is unlikely, if the masses of the modulus and saxion are close to
each other, the coupling of $aa\partial^{2}s_{r}$ induced by the kinetic
mixing can be much larger than $\kappa$. Such enhancement can make the
coupling of $s_{r}$ to NGBs suppressed not by $M_{Pl}$ but by $f$ around the
maximal mixing.

$\Delta{\cal L}|_{3}\ni-\frac{\kappa}{2}aa\partial^{2}s_{r}.$ (16)

The above shows that the modulus gets axion-philic, and decays into NGBs at
the rate

$\Gamma(s_{r}\to
aa)=\frac{\kappa^{2}}{32\pi}\frac{m^{3}_{s_{r}}}{M^{2}_{Pl}},$ (17)

where the masses of final states have been neglected. It should be noted that
non-suppression of the axionic decay mode $s_{r}\to aa$, which holds as long
as the modulus is lighter than or similar in mass to the scalar partner of the
NGB, is a model-independent generic feature of string moduli.

The non-derivative interactions involving $F$ and $F^{S}$ in (4) also mediate
the decay, which is however suppressed by the mass of final states or further.
This can be understood as follows. The matter $F$-component is expanded in
powers of $\phi$ as

$F=\langle
F\rangle+\langle\partial_{\phi}F\rangle\phi+\langle\partial_{\phi^{*}}F\rangle\phi^{*}+\cdots,$
(18)

while the scalar mass, $m_{\phi}$, is determined by

$m^{2}_{\phi}=\left\langle\frac{\partial^{2}V}{\partial\phi\partial\phi^{*}}\right\rangle=|\langle\partial_{\phi}F\rangle|^{2}+|\langle\partial_{\phi^{*}}F\rangle|^{2}+\cdots,$
(19)

where $V$ is the full scalar potential, and the ellipsis denotes possible
extra contributions. Hence, the coefficients of the linear terms in (18) are
bounded by

$|\langle\partial_{\phi}F\rangle|,\,|\langle\partial_{\phi^{*}}F\rangle|\lesssim
m_{\phi},$ (20)

barring cancellation among various contributions to the scalar mass. Expanded
in terms of $s$ and $s^{*}$, the modulus $F$-component exhibits the similar
feature.

Let us move on to stringy axions that appear in string compactification such
that a modulus, say $T$, is stabilized while preserving the associated shift
symmetry. The interactions responsible for the modulus decay come from the
Kähler potential

$(\delta S+\delta S^{*})(\delta T+\delta T^{*})^{2},$ (21)

where $\delta S$ and $\delta T$ are moduli fluctuations around the vacuum. The
$D$-term of the holomorphic cubic term, $\delta S\delta T\delta T$, includes
derivative interactions, which are however summed to be a total derivative.
The relevant derivative interactions from the above Kähler potential term read

$\Delta{\cal
L}|_{3}=\sqrt{2}\hat{\kappa}\left(\frac{1}{2}tt\partial^{2}s^{*}+st\partial^{2}t^{*}\right)+{\rm
h.c.},$ (22)

up to total derivatives. The coupling constant is determined by

$\hat{\kappa}=-\frac{1}{\sqrt{2}}\left\langle(\partial^{2}_{S}K_{0})^{-1/2}(\partial^{2}_{T}K_{0})^{-1}\partial_{S}\partial^{2}_{T}K_{0}\right\rangle,$
(23)

which is of the order unity in general. Here $t$ is the scalar component of
$\delta T$ and is decomposed as

$t=\frac{1}{\sqrt{2}}(\tau+i\varphi),$ (24)

for real scalars $\tau$ and $\varphi$. The stringy axion, $\varphi$, remains
massless until one adds non-perturbative effects breaking the shift symmetry
explicitly. It is easy to see that the derivative interactions are written in
terms of the real scalar fields as

$\displaystyle\hat{\kappa}^{-1}\Delta{\cal L}|_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\tau\tau\partial^{2}s_{r}-\frac{1}{2}\varphi\varphi\partial^{2}s_{r}+\varphi\tau\partial^{2}s_{i}$
(25) $\displaystyle+\,(\tau s_{r}-\varphi s_{i})\partial^{2}\tau+(\tau
s_{i}+\varphi s_{r})\partial^{2}\varphi.$

The first three interactions are expected to give a significant contribution
to the decays of $s_{r}$ and $s_{i}$, while the decay via the others is mass
suppressed.

To estimate correctly how fast the modulus decays via the derivative
interactions, we need to remove kinetic mixing induced by the Kähler potential
$(\delta S+\delta S^{*})(\delta T+\delta T^{*})$. Combined with the derivative
interactions from $(\delta S+\delta S^{*})^{3}$, $(\delta S+\delta
S^{*})^{2}(\delta T+\delta T^{*})$ and $(\delta T+\delta T^{*})^{3}$, the
kinetic mixing modifies the couplings of the moduli interactions involving
$\partial^{2}s_{r}$ and $\partial^{2}s_{i}$. However, their values in the
canonical mass basis generally remain the same order of magnitude because the
kinetic mixing is independent of $\hat{\kappa}$, differently from the NGB case
where the interactions (8) and (9) are all originated from the same single
Kähler potential term. The modulus decay into stringy axions is therefore
given by

$\Gamma(s_{r}\to\varphi\varphi)=\frac{\kappa^{\prime
2}}{32\pi}\frac{m^{3}_{s_{r}}}{M^{2}_{Pl}},$ (26)

where $\kappa^{\prime}$, which depends on $\hat{\kappa}$ and kinetic mixing
between the moduli, is generally of the order unity. Note that this non-
suppression of modulus decay into stringy axions is a quite general feature of
string compactification, and holds regardless of which one of $s_{r}$ and
$\tau$ is heavier.

The importance of modulus decay to stringy axions has been noticed for the
case where $S$ and $T$ are identical, in the the LARGE volume scenario
Cicoli:2012aq ; Higaki:2012ar , and later in the generalized setup
Higaki:2013lra . It is worth notifying that, if string compactification
involves a very light stringy axion, not only its scalar partner but also any
other moduli decay into stringy axions with a sizable branching fraction
insensitively to the details of moduli stabilzation. For instance, a massless
stringy axion arises in the generalized KKLT with multiple Kähler moduli if a
modulus, which does not appear in the superpotential, is stabilized by Kähler
mixing with the others Choi:2006za . In such a case, the decay into stringy
axions is sizable for all the real components of moduli.

###  II.2 Model-dependent decay modes

In this subsection we explore model-dependent interactions considerably
contributing to the modulus decay, and discuss how to make the modulus mainly
decay into the SM sector. Let us begin with a modulus coupling to the
holomorphic bilinear operator of matter superfields

$(\delta S+\delta S^{*})(\Phi_{1}\Phi_{2}+{\rm h.c.}),$ (27)

in the Kähler potential. It is the derivative interactions including only the
scalar fields that can enhance the modulus decay. Whereas the holomorphic
cubic term, $\delta S\Phi_{1}\Phi_{2}$, is irrelevant because its $D$-term
only gives a total derivative, the non-holomorphic cubic term gives

$\delta S^{*}\Phi_{1}\Phi_{2}|_{D}=-\phi_{1}\phi_{2}\partial^{2}s^{*},$ (28)

up to total derivatives. It is obvious that the decay of the modulus via the
interaction (28) is not suppressed by the mass of final states. Hence, an
interaction of the type (27) can enhance the modulus decay into the SM sector
if the involved matter fields are SM-charged.

More concretely, one can consider a modulus coupling to $QQ^{c}$ in the Kähler
potential with a coupling constant of the order unity or larger. Here the
matter superfields $Q+Q^{c}$ are vector-like under the SM gauge groups, and
are lighter than the modulus. A natural candidate of such operators is the
Higgs bilinear $H_{u}H_{d}$ Cicoli:2012aq ; Higaki:2012ar ; Higaki:2013lra ,
where $H_{u}$ and $H_{d}$ are respectively the up- and down-type Higgs
doublet. The modulus coupling to $H_{u}H_{d}$ in the Kähler potential

$\Delta K_{1}=\xi H_{u}H_{d}+{\rm h.c.}$ (29)

includes the term

$\langle\partial_{S}\xi\rangle(\delta S+\delta S^{*})H_{u}H_{d}+{\rm h.c.},$
(30)

where $\xi$ is a function of $S+S^{*}$. Meanwhile, it gives a contribution to
the Higgsino mass parameter $\mu$ and the Higgs mixing parameter $B$ as

$\Delta\mu=\langle\xi\rangle m_{3/2},\quad\Delta B=-m_{3/2},$ (31)

with $m_{3/2}$ being the gravitino mass, because it explicitly breaks the
super-Weyl symmetry. It is then difficult to avoid the Higgs $\mu/B\mu$
problem in the scenario where the gravitino is much heavier than the visible
sparticles, i.e. if anomaly mediation, which is a model-independent source of
soft supersymmetry breaking masses in supergravity, is sizable as is the case
in the KKLT moduli stabilization. To allow an order unity coupling of $\delta
S$ to $H_{u}H_{d}$ while suppressing the contributions to $\mu$ and $B$, one
can consider

$\Delta K_{2}=\xi\frac{X^{*}}{X}H_{u}H_{d}+{\rm h.c.},$ (32)

for a singlet $X$ that is radiatively stabilized by the loop potential
Nakamura:2008ey . As induced by the above super-Weyl invariant Kähler
potential term, both $\Delta\mu$ and $\Delta B$ are loop suppressed relative
to $m_{3/2}$ as required for the correct electroweak symmetry breaking, while
the modulus coupling to $H_{u}H_{d}$ is still given by
$\langle\partial_{S}\xi\rangle$ and is generally of the order unity.

Let us examine how the modulus decay is affected by the interaction
$(\ref{K-bilinear})$ when both $\phi_{1}$ and $\phi_{2}$ develop a VEV to
spontaneously break the global U$(1)$ symmetry, under which $\phi_{1}$ and
$\phi_{2}$ carry charge $1$ and $-1$, respectively, so that $\phi_{1}\phi_{2}$
becomes U$(1)$ invariant. The scalar fields are decomposed as

$\phi_{i}=\frac{f_{i}}{\sqrt{2}}\left(1+\frac{\sigma_{i}}{f_{i}}\right)e^{i\frac{a_{i}}{f_{i}}},$
(33)

for $i=1,2$. The massless NGB is then written

$a=\frac{f_{1}a_{1}-f_{2}a_{2}}{\sqrt{f^{2}_{1}+f^{2}_{2}}}.$ (34)

Through the interaction (28), the modulus couples to the two massive scalar
bosons composed of $\sigma_{i}$ and the other combination of $a_{i}$
proportional to $f_{2}a_{1}+f_{1}a_{2}$, but it is forbidden to couple to the
NGB due to the U$(1)$ symmetry. Nonetheless, in the presence of the coupling
(28), NGBs can be produced abundantly from the cascade decay of moduli because
the saxion generally decays mainly into NGBs via the $aa\partial^{2}\sigma$
interaction.

Another role of the interaction (28) is to induce kinetic mixing between the
modulus and $\sigma_{i}$ after spontaneous U$(1)$ breaking. Combined with the
interactions $a_{i}a_{i}\partial^{2}\sigma_{i}$ coming from the kinetic term
of $\phi_{i}$, such kinetic mixing generates the $aa\partial^{2}s_{r}$
interaction whose coupling is given by

$\frac{2f_{1}f_{2}}{f^{2}_{1}+f^{2}_{2}}$ (35)

times the coupling of (28), for the case where the modulus is much heavier
than the massive scalar bosons composed of $\sigma_{i}$. This implies that,
after electroweak symmetry breaking, the modulus coupling to $H_{u}H_{d}$ in
the Kähler potential induces the modulus decay into the CP even and odd
neutral Higgs bosons, and the charged Higgs bosons without mass suppression.
In addition, due to the kinetic mixing, it contributes to the modulus decay
into the associated NGBs, i.e. into the longitudinal modes of the weak gauge
bosons, according to the Goldstone boson equivalence theorem Cornwall:1974km ;
Vayonakis:1976vz

As a model-dependent modulus decay mode, one can also consider the decay into
the gauge sector, which is possible if the modulus appears in the gauge
kinetic function

$f_{a}=S+\Delta f_{a},$ (36)

with $\Delta f_{a}$ being a $S$-independent constant. It is straightforward to
see that the decay rate into gauge bosons reads Nakamura:2008ey

$\Gamma(s_{r}\to
gg)=\frac{N_{g}}{128\pi}\kappa^{2}_{g}\frac{m^{3}_{s_{r}}}{M^{2}_{Pl}},$ (37)

where $N_{g}=12$ counts the number of gauge bosons, and the order unity
constant $\kappa_{g}$ is determined by

$\kappa_{g}=\left\langle(\partial^{2}_{S}K_{0})^{-\frac{1}{2}}\frac{1}{{\rm
Re}f_{a}}\right\rangle.$ (38)

On the other hand, the decay rate into gauginos is written

$\Gamma(s_{r}\to\tilde{g}\tilde{g})=\beta^{2}_{r}\Gamma(s_{r}\to gg),$ (39)

with $\beta_{r}$ defined by

$\langle\partial_{s}F^{S}\rangle s+\langle\partial_{s^{*}}F^{S}\rangle
s^{*}\equiv\beta_{r}m_{s_{r}}s_{r}+\beta_{i}m_{s_{i}}s_{i},$ (40)

where $\beta_{r}$ and $\beta_{i}$ do not exceed order unity in size. Here the
decay rates have been evaluated neglecting the mass of decay products. The
decay rates of $s_{i}$ can be read off from those for $s_{r}$ by taking the
replacements, $m_{s_{r}}\to m_{s_{i}}$ and $\beta_{r}\to\beta_{i}$. Note that
the KKLT leads to $\beta_{r}\simeq\beta_{i}\simeq 1$ for $\Delta f_{a}=0$. In
order to enhance the decay into the visible sector further, one may consider
the case where the Higgsino mass parameter is generated from the Kähler
potential as (29) Cicoli:2012aq ; Higaki:2012ar  or as (32).

We close this section by mentioning that the modulus can decay to a gravitino
pair at a sizable rate if kinematically allowed. In particular, the decay is
dominated by the coupling to the helicity $\pm 1/2$ components, i.e. to the
Goldstino, which is proportional to the modulus $F$-term Endo:2006zj ;
Nakamura:2006uc ; Dine:2006ii . The decay rate into gravitinos reads

$\Gamma(s_{r}\to\tilde{G}\tilde{G})=\frac{\kappa^{2}_{3/2}}{288\pi}\frac{m^{3}_{s_{r}}}{M^{2}_{Pl}},$
(41)

in the limit $m_{3/2}\ll m_{s_{r}}$, with $m_{3/2}$ being the gravitino mass.
Here $\kappa_{3/2}$ is defined as follows

$\kappa_{3/2}\equiv\left\langle(\partial^{2}_{S}K_{0})^{\frac{1}{2}}F^{S}\right\rangle\frac{m_{s_{r}}}{m^{2}_{3/2}},$
(42)

whose size is below order unity because the modulus interactions are
gravitational. If it is of the order unity, as is typically the case in the
known examples of moduli stabilization, a large gravitino yield after moduli
decay can cause cosmological difficulties, dubbed as the moduli-induced
gravitino problem.

##  III Cosmological implications

Reheating takes place via moduli decay if the universe passes through a
modulus-dominated epoch as generally expected in string compactifications. To
catch the relevant features of such driven reheating, we simply consider a
single modulus $S$ under the assumption that its real and imaginary component
have a similar mass, $m_{s}$.333 For the case of multiple moduli, the decay
properties discussed in section II apply to all the moduli. If not much
heavier than the lightest one, heavy moduli can still be cosmologically
important. For instance, if there were two moduli, $S_{1}$ and $S_{2}$, with
masses $m_{1}$ and $m_{2}$ ($m_{1}<m_{2}$), respectively, the energy density
of the decay products of $S_{2}$ is diluted by the entropy released from late-
time decay of $S_{1}$. The dilution factor is given by the ratio between their
decay temperatures, and reads $\Delta\sim(m_{2}/m_{1})^{3/2}$, assuming that
the displacement of moduli at the onset of coherent oscillation is of the
order of the Planck scale.  In flux compactifications, a modulus unfixed by
flux can be stabilized by non-perturbative dynamics or Kähler mixing with
others. Its mass is then tied to the scale of supersymmetry breaking,
$m_{s}\sim m_{3/2}$, up to a factor of order $\ln(M_{Pl}/m_{3/2})$ Choi:2005ge
; Choi:2008hn . On the other hand, in the LARGE volume scenario, the mass of
the large volume modulus can be much smaller than $m_{3/2}$
Balasubramanian:2005zx .

To allow successful BBN after the modulus dominated phase, the modulus should
decay dominantly into the visible sector while forming the standard thermal
background with temperature $T_{s}$ constrained by deSalas:2015glj ;
Hasegawa:2019jsa

$T_{s}\simeq\left(\frac{90}{\pi^{2}g_{\ast}}\right)^{\frac{1}{4}}(\Gamma_{s}M_{Pl})^{\frac{1}{2}}\gtrsim
5\,{\rm MeV},$ (43)

where $g_{\ast}$ is the relativistic degrees of freedom at $T_{s}$, and
$\Gamma_{s}$ is the total decay width of the modulus. Here we have assumed an
instantaneous conversion of the modulus energy density to radiation. The
constraint on $T_{s}$ requires the modulus to be heavier than about $100$ TeV.
Modulus domination also constrains models for baryon and dark matter genesis
because a huge amount of entropy is released from moduli decay. If the genesis
occurs at the early stage of moduli domination or before, the dilution factor
naively reads, $\Delta\sim 10^{12}\,(m_{s}/{\rm PeV})^{-1}$, since the initial
displacement of the modulus after inflation is generally of the order of the
Planck scale.

A viable reheating requires the modulus to decay dominantly to the visible
sector. This can be achieved by adding model-dependent modulus couplings such
as a coupling to the bilinear operator $QQ^{c}$ or $H_{u}H_{d}$ in the Kähler
potential, or to the visible gauge sector as discussed in subsection II.2. In
the presence of such couplings, the decay rate into the visible sector reads

$\Gamma_{s\to{\rm SM}}=\frac{\kappa^{2}_{\rm
SM}}{32\pi}\frac{m^{3}_{s}}{M^{2}_{Pl}},$ (44)

where $\kappa_{\rm SM}$ is a constant of order unity. As discussed already,
the modulus generally decays into axions without mass suppression if string
compactification involves a light stringy axion, $\varphi$, or if there is a
(pseudo) NGB, $a$, whose scalar partner has a mass larger than or comparable
to the modulus mass. A theoretically well-motivated NGB is the QCD axion
Ringwald:2012hr . The branching ratio of the axion-philic modulus into axions
is written

${\rm Br}(s\to{\rm axions})=\frac{\kappa^{2}+\kappa^{\prime
2}}{\kappa^{2}_{\rm SM}},$ (45)

which naturally lies in the range of the order of $0.1$. If the axions are
cosmologically stable, the produced axions form dark radiation. The
contribution to the effective number of neutrino species is estimated by
Choi:1996vz

$\Delta N_{\rm eff}=\frac{43}{7}\frac{{\rm Br}(s\to{\rm axions})}{1-{\rm
Br}(s\to{\rm axions})}\left(\frac{g_{\ast}(T_{\nu,{\rm
dec}})}{g_{\ast}(T_{s})}\right)^{\frac{1}{3}},$ (46)

under the assumption that the branching ratio into the visible sector is given
by $1-{\rm Br}(s\to{\rm axions})$. Here $T_{\nu,{\rm dec}}$ is the neutrino
decoupling temperature. An observable amount of dark radiation is a natural
prediction of the axion-philic modulus insensitively to the detail of moduli
stabilization. Dark radiation can have interesting cosmological effects in the
early universe. For instance, if $\Delta N_{\rm eff}\gtrsim 0.5$, which is the
case for ${\rm Br}(s\to{\rm axions})$ larger than about $0.07$, it may relieve
the Hubble tension Riess:2019cxk ; DiValentino:2021izs .

It should be noted that the presence of a plenitude of light stringy axions,
an axiverse, has been suggested as evidence for the extra dimensions of string
theory Arvanitaki:2009fg . In the axiverse, the lightest modulus would
dominantly decay into these axions, producing too much dark radiation if the
universe undergoes a modulus-dominated phase. The modulus coupling to
$H_{u}H_{d}$ or to the gauge sector would not be sufficient to suppress ${\rm
Br}(s\to{\rm axions})$. To avoid the moduli-induced axion problem, i.e. to
enhance the modulus decay into the visible sector, one may thus rely on the
modulus coupling to $QQ^{c}$ in the Kähler potential, for instance, for a
number of $Q+Q^{c}$ that form complete SU$(5)$ multiplets. Here the
supersymmetric mass of $Q+Q^{c}$ in the superpotential is assumed to be
smaller than the modulus mass. We also note that, depending on their masses,
the number of $Q+Q^{c}$ is constrained by the perturbativity of the gauge
interactions up to the grand unification scale. The phenomenological and
cosmological details of such particles are out of the scope of this paper, and
hence will not be touched here.

The cosmological features discussed so far remain the same even if the saxion
has an initial amplitude of the order of the Planck scale, as long as it
decays much faster than the modulus, i.e. if

$m_{\sigma}\gtrsim 100\,{\rm GeV}\left(\frac{f}{10^{12}{\rm
GeV}}\right)^{\frac{2}{3}}\left(\frac{m_{s}}{\rm PeV}\right).$ (47)

However, the situation changes much if the saxion is thermally trapped at the
origin to drive thermal inflation Lyth:1995hj . It is also a plausible
scenario although we do not consider it in this paper.

Moduli stabilization with $m_{s}>2m_{3/2}$ generally confronts the moduli-
induced gravitino problem unless one adds $R$-parity violating interactions.
Interestingly, the gravitino problem can be relaxed if the fermionic partner
of the axion, the axino, is much lighter than the lightest observable
sparticle (LOSP) since otherwise LOSPs produced from gravitino decay would
overclose the universe under the $R$-parity conservation. Even for moduli
stabilization with $m_{s}<2m_{3/2}$, the axino is still an attractive particle
because, as a FIMP, it is a viable alternative to weakly interacting massive
particle (WIMP) cold dark matter.

In the moduli-dominated universe, axinos are non-thermally produced from the
decays of the modulus, gravitino, and LOSP. The axino abundance from LOSPs
decay is quite model-dependent because LOSPs can efficiently annihilate before
their decay depending on the LOSP nature and the value of $f$. Meanwhile,
independently of model details, the axino yield directly produced from the
gravitationally interacting moduli and gravitinos is estimated by

$Y_{\tilde{a}}|_{s\to\tilde{a}}=2B^{s}_{\tilde{a}}\frac{3}{4}\frac{T_{s}}{m_{s}},$
(48)

where $B^{s}_{\tilde{a}}$ is defined by

$B^{s}_{\tilde{a}}\equiv{\rm Br}(s\to\tilde{a}\tilde{a})+{\rm
Br}(s\to\tilde{G}\tilde{G}){\rm Br}(\tilde{G}\to\tilde{a}),$ (49)

with ${\rm Br}(s\to\tilde{a}\tilde{a})\sim(m_{\tilde{a}}/m_{s})^{2}$ as
follows from (4). Thus, in the case with $m_{s}<2m_{3/2}$, the axino should
have mass,

$m_{\tilde{a}}\lesssim 1\,{\rm TeV}\left(\frac{0.1{\rm
GeV}}{T_{s}}\right)^{\frac{1}{3}}\left(\frac{m_{s}}{\rm PeV}\right).$ (50)

not to exceed the observed dark matter density. On the other hand, if
$m_{s}>2m_{3/2}$, the gravitino produced moduli uniformly decays to the axino
and all lighter visible sparticles, and should be heavier than a few tens of
TeV in order to decay before BBN Kawasaki:2017bqm . The branching fraction of
the modulus decay to gravitinos is generally sizable, indicating that the
axino should be much lighter than the above bound. The axino mass and the LOSP
nature may further be constrained by the axino free-streaming scale at the
time of matter-radiation equality, which should be less than the order of Mpc
to be consistent with the observations of Lyman-$\alpha$ forest Viel:2005qj .
Axinos are mostly produced from the decay of LOSPs as long as the number of
the visible sparticles lighter than the gravitino is much larger than unity.
Then, the free-streaming constraint is relaxed if the LOSP scattering rate is
larger than its decay rate so that LOSPs become non-relativistic before their
decay. See refs. Nakamura:2008ey ; Bae:2021rmg  for some examples in the KKLT
scenario. Finally, we note that axinos are also produced via freeze-in
processes, which can be effective depending of the value of $f$ if $T_{s}$ is
not much lower than the freeze-out temperature of the LOSP Bae:2021rmg .

##  IV Conclusions

In this work we have showed that string moduli have axion-philic nature.
Induced by the model-insensitive derivative interactions, the decay of a
modulus into stringy axions occurs without suppression by the mass of decay
products. Interestingly, we found that this feature holds not only for the
scalar partner of the stringy axion but also for any other moduli. The decay
into (pseudo) NGBs also avoids mass suppression if the modulus is not much
heavier than the saxion. The axion-philic nature makes string moduli a natural
source of an observable amount of dark radiation in the universe that passes
through a modulus-dominated era as is generally expected in scenarios of
string compactification. Proper reheating would then require moduli to have
additional couplings to enhance their decay into the SM sector so that hot
thermal bath composed of mostly SM particles is formed. For instance, the
modulus can couple to the gauge sector if it appears in the gauge kinetic
function. Another way is to add a modulus coupling to the holomorphic bilinear
operator of matter superfields in the Kähler potential. It is also interesting
to note that the fermionic superpartner of a NGB as well as coherently
oscillating NGBs can be a viable alternative to WIMP dark matter.

Acknowledgments

The authors thank to Fuminobu Takahashi for helpful discussion. This work was
supported by the National Research Foundation of Korea (NRF) grant funded by
the Korea government: 2018R1C1B6006061, 2021R1A4A5031460 (K.S.J.),
2017R1D1A1B06035959 (W.I.P.), and also by Research Base Construction Fund
Support Program funded by Jeonbuk National University in 2021 (W.I.P.).

## References

  * (1) G. D. Coughlan, W. Fischler, E. W. Kolb, S. Raby and G. G. Ross, Phys. Lett. B 131 (1983), 59-64.
  * (2) J. R. Ellis, D. V. Nanopoulos and M. Quiros, Phys. Lett. B 174 (1986), 176-182.
  * (3) A. S. Goncharov, A. D. Linde and M. I. Vysotsky, Phys. Lett. B 147 (1984), 279-283.
  * (4) M. Endo, K. Hamaguchi and F. Takahashi, Phys. Rev. Lett. 96 (2006), 211301 [arXiv:hep-ph/0602061 [hep-ph]].
  * (5) S. Nakamura and M. Yamaguchi, Phys. Lett. B 638 (2006), 389-395 [arXiv:hep-ph/0602081 [hep-ph]].
  * (6) M. Dine, R. Kitano, A. Morisse and Y. Shirman, Phys. Rev. D 73, 123518 (2006) [hep-ph/0604140].
  * (7) M. Cicoli, J. P. Conlon and F. Quevedo, Phys. Rev. D 87 (2013) no.4, 043520 [arXiv:1208.3562 [hep-ph]].
  * (8) T. Higaki and F. Takahashi, JHEP 11 (2012), 125 [arXiv:1208.3563 [hep-ph]].
  * (9) T. Higaki, K. Nakayama and F. Takahashi, JHEP 07, 005 (2013) [arXiv:1304.7987 [hep-ph]].
  * (10) S. Kachru, R. Kallosh, A. D. Linde and S. P. Trivedi, Phys. Rev. D 68 (2003), 046005 [arXiv:hep-th/0301240 [hep-th]].
  * (11) V. Balasubramanian, P. Berglund, J. P. Conlon and F. Quevedo, JHEP 03 (2005), 007 [arXiv:hep-th/0502058 [hep-th]].
  * (12) S. Nakamura, K. i. Okumura and M. Yamaguchi, Phys. Rev. D 77 (2008), 115027 [arXiv:0803.3725 [hep-ph]].
  * (13) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977); Phys. Rev. D 16, 1791 (1977).
  * (14) K. Choi and K. S. Jeong, JHEP 01 (2007), 103 [arXiv:hep-th/0611279 [hep-th]].
  * (15) J. M. Cornwall, D. N. Levin and G. Tiktopoulos, Phys. Rev. D 10 (1974), 1145 [erratum: Phys. Rev. D 11 (1975), 972].
  * (16) C. E. Vayonakis, Lett. Nuovo Cim. 17 (1976), 383.
  * (17) K. Choi, A. Falkowski, H. P. Nilles and M. Olechowski, Nucl. Phys. B 718 (2005), 113-133 [arXiv:hep-th/0503216 [hep-th]].
  * (18) K. Choi, K. S. Jeong and K. I. Okumura, JHEP 07 (2008), 047 [arXiv:0804.4283 [hep-ph]].
  * (19) P. F. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor and O. Pisanti, Phys. Rev. D 92 (2015) no.12, 123534 [arXiv:1511.00672 [astro-ph.CO]].
  * (20) T. Hasegawa, N. Hiroshima, K. Kohri, R. S. L. Hansen, T. Tram and S. Hannestad, JCAP 12 (2019), 012 [arXiv:1908.10189 [hep-ph]].
  * (21) For recent reviews, see, for instance, A. Ringwald, Phys. Dark Univ. 1 (2012), 116-135 [arXiv:1210.5081 [hep-ph]]; M. Kawasaki and K. Nakayama, Ann. Rev. Nucl. Part. Sci. 63 (2013), 69-95 [arXiv:1301.1123 [hep-ph]].
  * (22) K. Choi, E. J. Chun and J. E. Kim, Phys. Lett. B 403 (1997), 209-217 [arXiv:hep-ph/9608222 [hep-ph]].
  * (23) A. G. Riess, S. Casertano, W. Yuan, L. M. Macri and D. Scolnic, Astrophys. J. 876 (2019) no.1, 85 [arXiv:1903.07603 [astro-ph.CO]].
  * (24) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess and J. Silk, Class. Quant. Grav. 38 (2021) no.15, 153001 [arXiv:2103.01183 [astro-ph.CO]].
  * (25) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper and J. March-Russell, Phys. Rev. D 81 (2010), 123530 [arXiv:0905.4720 [hep-th]].
  * (26) D. H. Lyth and E. D. Stewart, Phys. Rev. Lett. 75 (1995), 201-204 [arXiv:hep-ph/9502417 [hep-ph]]; Phys. Rev. D 53 (1996), 1784-1798 [arXiv:hep-ph/9510204 [hep-ph]].
  * (27) M. Kawasaki, K. Kohri, T. Moroi and Y. Takaesu, Phys. Rev. D 97 (2018) no.2, 023502 [arXiv:1709.01211 [hep-ph]].
  * (28) M. Viel, J. Lesgourgues, M. G. Haehnelt, S. Matarrese and A. Riotto, Phys. Rev. D 71 (2005), 063534 [arXiv:astro-ph/0501562 [astro-ph]].
  * (29) K. J. Bae and K. S. Jeong, Phys. Rev. D 104 (2021) no.1, 015013 [arXiv:2105.08236 [hep-ph]].


11institutetext: Université Paris-Saclay, CNRS, Institut d’Astrophysique
Spatiale, Bâtiment 121, 91405 Orsay, France 22institutetext: Department of
Physics and Astronomy, University of British Columbia, Vancouver, BC V6T 1Z1,
Canada 33institutetext: Canadian Institute for Advanced Research, 180 Dundas
St W, Toronto, ON M5G 1Z8, Canada 44institutetext: Canada Research Chair in
Observational Cosmology 55institutetext: Astrophysics Research Institute,
Liverpool John Moores University, Liverpool, L3 5RF, United Kingdom
66institutetext: Astrophysics and Cosmology Research Unit, School of Chemistry
and Physics University of KwaZulu-Natal, Durban, South Africa 77institutetext:
NAOC-UKZN Computational Astrophysics Center (NUCAC), University of KwaZulu-
Natal, Durban, 4000 88institutetext: Institute for Astronomy, University of
Edinburgh, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, UK
99institutetext: Department of Physics and Astronomy, University College
London, Gower Street, London, WC1E 6BT, UK 1010institutetext: Instituto de
Fisica, Universidade Federal do Rio de Janeiro, 21941-972, Rio de Janeiro, RJ,
Brazil

# Constraining AGN feedback model with SZ profile

Hideki Tanimura 11<EMAIL_ADDRESS>Gary Hinshaw 223344 Ian G.
McCarthy 55 Ludovic Van Waerbeke 22 Nabila Aghanim Yin-Zhe Ma 116677 Alexander
Mead 88 Tilman Tröster 88 Alireza Hojjati 22 Bruno Moraes 991010

###### Abstract

Relativistic jets from AGN have a wide range of impacts on galaxy groups and
clusters and are key for understanding their formation and physical
properties. However, this non-gravitational process is not well understood.
Galaxy groups with shallow gravitational potentials are ideal laboratories to
study and constrain the AGN feedback model. We studied hot gas in $\sim$66,000
SDSS LRG halos with an average halo mass of $3\times 10^{13}$ M⊙ using the
Planck tSZ map. We have detected their average tSZ radial profile at
$\sim$17$\sigma$ and compared it with the cosmo-OWLS cosmological
hydrodynamical simulations with different AGN feedback models. The best
agreement has been obtained for the AGN 8.0 model in the simulations. We have
also compared our measured tSZ profile with the prediction from the universal
pressure profile assuming the self-similar relation and found them consistent
if the model accounts for the clustering of neighboring haloes via a two-halo
term.

##  1 Introduction

A power-law relation between the electron (or gas) pressure and mass,
$P-M^{2/3}$, is valid under the assumption that the galaxy-formation process
is dominated by gravity, called self-similar relation. Any deviation from this
relation points to the presence of more complex processes such as baryonic
feedback effects. X-ray and SZ observations demonstrate that massive galaxy
clusters almost follow the self-similar scaling relation Planck2013b . On the
other hand, numerous X-ray observations find a deficit of baryonic gas inside
low-mass halos compared to the cosmic average baryon fraction Gastaldello2007
; Pratt2009 ; Sun2009 ; Gonzalez2013 , implying the deviation from the self-
similar scaling relation. The gas deficit may be caused by AGN feedback: the
response of the intracluster medium (ICM) to the relativistic plasma ejected
by the active galactic nuclei (AGN) Gitti2012 . The AGN feedback is key for
understanding the evolution of galaxies and the formation of large-scale
structures; however, not well understood.

##  2 Stacking $y$ map centered on LRGs

We performed a stacking analysis using the Planck $y$ map Planck2016-XXII  for
Sloan Digital Sky Survey (SDSS) data release 7 (DR7) luminous red galaxies
(LRGs) at $0.16<z<0.47$ Blanton2005 , which are considered to be mostly
central galaxies in dark matter halos Reid2009 . We constructed an average $y$
profile centered on the 66,479 LRGs with the stellar mass111The stellar masses
of the LRGs are provided in the New York University Value-Added catalogue
(NYU-VAGC) at http://sdss.physics.nyu.edu/vagc/ of $10^{11.2}\leq
M_{*}/M_{\odot}\leq 10^{11.7}$, roughly corresponding to halo masses of
$10^{13}\leq M_{500}/M_{\odot}\leq 10^{14}$. We detected a significant signal
at $\sim 17\sigma$ out to $\sim$30 arcmins well beyond the extent of the 10
arcmin beam of the Planck $y$ map in Fig. 1.

Figure 1:  Left: Average $y$ map centered on 66,479 LRGs. Right: Average $y$
profile around the LRGs. The 1$\sigma$ statistical uncertainty is represented
via the width of the blue line. The $\mathrm{FWHM}=10$ arcmin Gaussian beam of
the Planck $y$ map is shown as a black dashed line for comparison, the peak of
which is normalized to the center of the LRGs’ $y$ profile.

##  3 Comparison with the cosmo-OWLS hydrodynamic simulations

We compared the average $y$ profile around LRGs with predictions from the
cosmo-OWLS cosmological hydrodynamical simulations Brun2014 ; Daalen2014 ;
McCarthy2014 . The simulations are designed with cluster cosmology and large-
scale structure surveys in mind. They consist of box-periodic hydrodynamical
simulations and have volumes of 400$h^{-1}\rm{Mpc}^{3}$ with $1024^{3}$
baryonic and dark matter particles. They simulated five different models
(Table. 1), which are compared with our measured $y$ profile in Fig. 2.

In the comparison, a clear difference between the data and REF model can be
seen. In general, energy released from the center of a halo heats cluster gas,
this in turn prevents cooling and thus the star formation around the central
region. Therefore, if we consider halos of the same total mass, the stellar
mass of the central galaxy is decreased as the power of the central AGN is
increased. Since we select central galaxies based on stellar mass, lower-mass
halos are selected in the REF model compared to the models that include AGN
feedback. This is apparent as the lower central peak value of the simulated
$y$ profiles in the REF model compared to the AGN models. We also see a
visible trend in that the higher the power of AGN feedback, the lower the peak
of $y$ profile. This is due to the fact that the AGN feedback ejects gas from
the center of halos outward and the overall gas density is lowered.

In summary, this comparison agrees best with simulations that include AGN
feedback (AGN 8.0) but not with simulations that do not include it (REF) or
with simulations with very violent AGN feedback (AGN 8.7). The reduced
$\chi^{2}$ values are 1.0 (AGN 8.0), 1.3 (AGN 8.5), 9.9 (AGN 8.7), and 75.9
(REF), respectively. This result is consistent with other studies showing that
the AGN 8.0 model reproduces a variety of observed gas features in optical and
X-ray data (e.g., Brun2014 ).

Table 1: The baryon feedback models in the cosmo-OWLS simulations. Simulation | Cooling | Star formation | SN feedback | AGN feedback | $\Delta T_{\rm heat}$
---|---|---|---|---|---
NOCOOL | No | No | No | No | …
REF | Yes | Yes | Yes | No | …
AGN 8.0 | Yes | Yes | Yes | Yes |  $10^{8.0}$ K
AGN 8.5 | Yes | Yes | Yes | Yes |  $10^{8.5}$ K
AGN 8.7 | Yes | Yes | Yes | Yes |  $10^{8.7}$ K

Figure 2: The average $y$ profile around LRGs (blue) is compared to the $y$
profiles of the simulated central galaxies (gray) in different AGN feedback
models respectively. In each case we have matched the stellar mass and
redshift distributions from the simulations to be the same as those in the
data. Top left: AGN 8.0 model, Top right: AGN 8.5 model, Bottom left: AGN 8.7
model, and Bottom right: REF model.

##  4 Halo model prediction with the UPP

The average $y$ profile around the LRGs is also compared with a prediction
using the halo model with the halo mass function and halo bias Tinker2010  and
a universal pressure profile (UPP) Planck2013-V  in Fig. 3. The stellar masses
of the LRGs are converted to halo masses using the stellar-halo mass relation
from either Coupon2015  (C15-SHM) or Wang2016  (W16-SHM), estimated using
gravitational lensing measurements. Then, the UPP is scaled to different
masses of halos, assuming a self-similar relation. The predicted $y$ profile
is consistent with the data, but only if we account for the two-halo
clustering term in the model and if we assume the stellar-halo mass relation
from C15-SHM or W16-SHM. This consistency may imply that the UPP, estimated
for massive galaxy clusters in the mass range of $10^{14.4}-10^{15.3}$ M⊙, can
be applied even in low-mass halos down to $\sim 10^{13.5}$ M⊙ in a self-
similar manner.

##  5 Conclusion

According to the study with the cosmo-OWLS simulations Brun2015 , the AGN 8.0
simulation predicts more extended pressure profiles around low-mass halos than
the UPP. However, the AGN 8.0 model also reproduces our observed $y$ profile,
supporting a self-similar assumption. This apparent inconsistency is explained
by the fact that the deviations from the self-similar relation in the AGN 8.0
simulation are mainly confined to halo masses below $M_{500}\sim 10^{13.5}$
M⊙, which roughly corresponds to the average mass of our sample. Furthermore,
the impact of coarse angular resolution of the Planck $y$ map is not
negligible in our analysis: the UPP and AGN 8.0 pressure distributions only
differ significantly on scales of $r<r_{500}$ in the simulations, in which
they are compared without the beam. Thus, data from higher-resolution tSZ maps
(such as those from ACT or SPT with FWHM of order an arcminute) would be
necessary in this regard.

Figure 3: The average $y$ profile around 66,479 LRGs (blue) is compared to the
predictions using a halo model with the halo mass function and halo bias
Tinker2010  and UPP Planck2013-V . The halo masses of the LRGs are estimated
using either the SHM relation of Coupon2015  (C15-SHM: magenta) and Wang2016
(W16-SHM: yellow). The one-halo (dashed line in magenta) and two-halo (dash-
dotted line in magenta) terms are shown separately for the model prediction
using the C15-SHM. The $y$ profile of the simulated central galaxies in the
AGN 8.0 simulation is shown in grey. To show an impact of beam, the average
angular size of the LRGs ($\bar{\theta}_{500}\sim 1.6$ arcmin) is shown in
vertical black dashed line.

## Acknowledgement

This research has been supported by the funding for the ByoPiC project from
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme grant agreement ERC-2015-AdG 695561.

## References

  * (1) Planck Collaboration, A&A, 557, A52 (2013)
  * (2) Gastaldello, F., et al., ApJ, 669, 158 (2007)
  * (3) Pratt, G. W., et al., A&A, 498, 361 (2009)
  * (4) Sun, M., et al., ApJ, 693, 1142 (2009)
  * (5) Gonzalez, A. H., et al., The Astrophysical Journal, 778, 14 (2013)
  * (6) Gitti, Myriam, et al., Advances in Astronomy, 2012, 950641 (2012)
  * (7) Planck Collaboration, A&A, 594, A22 (2016)
  * (8) Blanton, Michael R., et al., AJ, 129, 2562-2578 (2005)
  * (9) Reid, B. A., et al., AJ, 698, 143-154 (2009)
  * (10) Le Brun, Amandine M. C., et al., MNRAS, 440, 1270-1290 (2014)
  * (11) van Daalen, M. P., et al., MNRAS, 441, 2997-3010 (2014)
  * (12) McCarthy, I. G., et al., MNRAS, 440, 3645-3657 (2014)
  * (13) Tinker, Jeremy L., et al., ApJ, 724, 878-886 (2010)
  * (14) Planck Collaboration, A&A, 550, A131 (2013)
  * (15) Coupon, J., et al., MNRAS, 449, 1352-1379 (2015)
  * (16) Wang, W., et al., MNRAS, 456, 2301-2320 (2016)
  * (17) Le Brun, Amandine M. C., et al., MNRAS, 451, 3868-3881 (2015)


# Many Body Density of States of a system of non interacting spinless fermions

Rémi Lefèvre  Krissia Zawadzki  Grégoire Ithier  Department of Physics, Royal
Holloway University of London<EMAIL_ADDRESS>
###### Abstract

The modeling of out-of-equilibrium many-body systems requires to go beyond the
low-energy physics and local densities of states. Many-body localization,
presence or lack of thermalization and quantum chaos are examples of phenomena
in which states at different energy scales, including the highly excited ones,
contribute to the dynamics and therefore affect the system’s properties.
Quantifying these contributions requires the many-body density of states
(MBDoS), a function whose calculation becomes challenging even for non-
interacting identical quantum particles due to the difficulty in enumerating
states while enforcing the exchange symmetry. In the present work, we
introduce a new approach to evaluate the MBDoS in the case of systems that can
be mapped into free fermions. The starting point of our method is the
principal component analysis of the filling matrix $F$ describing how $N$
fermions can be configured into $L$ single-particle energy levels. We show
that the many body spectrum can be expanded as a weighted sum of spectra given
by the principal components of the filling matrix. The weighting coefficients
only involve renormalized energies obtained from the single body spectrum. We
illustrate our method in two classes of problems that are mapped into spinless
fermions: (i) non-interacting electrons in a homogeneous tight-binding model
in 1D and 2D, and (ii) interacting spins in a chain under a transverse field.

##  I Introduction

The concept of Density of States (DoS) is at the heart of statistical physics
where it defines partition function and temperature. In nuclear physics,
quantifying the level density is necessary to describe nuclear reactions
involving excited states Zelevinsky and Horoi (2019). No less importantly, it
is one of the most appealing quantity in condensed matter physics, where one
is interested in investigating how electrons and holes populate energy bands
to give rise to material’s properties Ashcroft  _et al._ (1976). In all these
fields, the density of states is crucial to characterize a multi-particle
system and determine which states are accessible at energy scales of interest.

For a long time, the success of mean field theories and the quasi-particle
picture to describe a highly degenerate Fermi liquid has promoted the Single-
Body Density of States (SBDoS) and its local counterpart to the focus of
investigations of electronic systems. In the presence of interactions, efforts
have been concentrated in the physics at low temperatures, so that the SBDoS
around the Fermi level suffices to obtain most of their properties.
Nonetheless, the quest of calculating a Many-Body Density of States (MBDoS)
has became arguably necessary in the context of isolated quantum systems
undergoing an out-of-equilibrium dynamics Gogolin and Eisert (2016). There,
contributions from different parts of the spectrum prevents one from relying
on a description based solely on the low-lying states and on the SBDoS.
Quantifying these contributions is crucial to shed light on phenomena such as
quantum chaos Santos and Rigol (2010), thermalization and its lack Borgonovi
_et al._ (2016); Polkovnikov  _et al._ (2011); Nandkishore and Huse (2015),
many-body localization Abanin  _et al._ (2019) and more generally
unconventional stationary states Ithier and Benaych-Georges (2017); Ithier
_et al._ (2017).

In this respect, a MBDoS provides useful information, as it allows for
quantifying how interactions between individual constituents lead to a complex
many-body dynamics, with coexisting single-particle and collective effects.
Even non-interacting systems pose a challenge due to the combinatorial nature
of the problem of how single-body levels can be populated to define a
distribution of many-body energy levels. For both non-interacting and
interacting systems, it is possible to retrieve numerically the full spectrum
of energies and eigenstates only for system sizes that do not exceed a dozen
of particles. Symmetries can aid this computation by allowing one to split the
total Hilbert space into blocks, associated with conserved quantum numbers.
This idea has served as the basis of exact diagonalization Weiße and Fehske
(2008) and it is also implemented in well established numerical methods, as
for instance the Kernel Polynomial Method (KPM) Weiße  _et al._ (2006). On the
analytical side, methods to calculate the MBDoS started in nuclear physics
with a Fermi gas approximation calculation derived by Bethe Bethe (1936),
which inspired approaches such as the constant temperature or the continuum
shell model (see e.g. the reviews in Zelevinsky and Horoi (2019); Volya and
Zelevinsky (2006)). More involved methods using exact combinatorial counting
Hillman and Grover (1969); Berger and Martinot (1974), recursive relations
Jacquemin and Kataria (1986), or saddle approximations Bohr and Mottelson
(1998) provided some approximate results.

The calculation of the MBDoS of systems of non interacting spinless fermions
is the problem at the focus of the present paper. We propose a new approach to
calculate the exact MBDoS based on the symmetries of a rectangular filling
matrix describing how $N$ particles can be combined into $L$ single-particle
energy levels to generate the many-body states. The starting point of our
method relies on the singular value decomposition of this matrix, which allows
to expand the many-body spectrum as a weighted sum of ’principal’ spectra.
These spectra depend only on the number of single body levels and the number
of particles. The weighting factors involve a discrete Fourier transformation
of the single body energies, providing renormalized energies. Our approach is
illustrated in the case of spinless non-interacting fermions, and applied to
the tight-binding model in one and two dimensions, and to the transverse Ising
field chain.

We organize the paper as follows. In Sec. II, we introduce the idea of
principal analysis decomposition of a filling matrix and explain how it allows
to access a non interacting Many Body spectrum. In Sec. III, we develop the
method in the case of spinless fermions, discuss how to explore symmetries of
the problem of calculating the MBDoS. Applications to tight-binding and Ising
chains are discussed in Sec. IV. Finally, our main findings are summarized in
Sec. V.

##  II Principal Components Approach

We start by considering a system for which we know the single-body spectrum
given by the energies $\epsilon_{k}$. Our goal is to construct the many-body
spectrum, which has energies $E_{p}=\sum_{k=1}^{L}\epsilon_{k}f_{p}^{k}$,
where $f_{p}^{k}$ is the filling factor or occupation number of the $k^{th}$
single body level for the $p$-th many-body state. Our first step is to rewrite
all those energies in matrix form by collecting all single-body energies
$\epsilon_{k}$ into a column vector $\epsilon$, all many body energies $E_{p}$
into a column vector $E$, and finally construct a rectangular matrix $F$, the
”filling” matrix, of general term $f_{p}^{k}$. We can then write the relation
between many body and single body energies in matrix form :

$E=F\cdot\epsilon$ (1)

where $\cdot$ denotes the matricial product.

The singular value decomposition (SVD) of the filling matrix
$F=U\cdot\Sigma\cdot V^{\dagger}$ with $U,V$ square unitary matrices and
$\Sigma$ rectangular diagonal, provides its so called ”principal component”
expansion :

$F=\sum_{\ell=0}^{L-1}\sigma_{\ell}\;U_{\ell}\cdot V_{\ell}^{\dagger}$ (2)

where $U_{\ell}$ is the $l^{th}$ right-singular vector (i.e. the $l^{th}$
column of $U$), $\sigma_{\ell}$ are the singular values (i.e. the main
diagonal of $\Sigma$), and $V_{\ell}$ is the $l^{th}$ left-singular vector
(i.e. the $l^{th}$ column of $V$). To ease notations, in the following we will
incorporate the singular value $\sigma_{l}$ into the definition of $U_{l}$ so
that the right-singular vectors are orthogonal but not normalized anymore.
Note that with these notations, $U_{\ell}$ and $V_{\ell}$ are simply related
by $U_{\ell}=F\cdot V_{\ell}$. By combining Eqs. 1 and 2, the vector $E$ of
many-body energies can be rewritten as :

$E=\sum_{\ell=0}^{L-1}U_{\ell}(V_{\ell}^{\dagger}\cdot\epsilon)=\sum_{\ell=0}^{L-1}\tilde{\epsilon}_{\ell}\;U_{\ell}$
(3)

where we defined the scalars
$\tilde{\epsilon_{\ell}}=V_{\ell}^{\dagger}\cdot\epsilon$.

The many body spectrum now appears as a sum of ”principal” spectral components
given by the $U_{\ell}$ vectors. These spectral components are weighted by the
new effective energies $\tilde{\epsilon_{\ell}}$ of a renormalized single-body
spectrum. In particular, the form of $V_{\ell}$ defines which of these
energies contribute to the many-body spectrum. As we will discuss later, a
very convenient choice is given by Fourier modes obtained from the analytical
solution for the eigenvectors of the circulant matrix $F^{\dagger}\cdot F$.
Depending on the band structure of the system, it is also possible that some
renormalized energies vanish which allows the associated $\ell$ modes to be
discarded.

The problem of computing the many body spectrum is now split in two parts. The
first part is to compute the renormalized single body spectrum only depending
on the right singular part of the SVD (ie. the $V_{\ell}$) which can be
obtained analytically or for a small computational cost scaling like $L$. The
second part depends on the left singular part of the SVD (ie. the $U_{\ell}$
incorporating the $\sigma_{l}$) which only contain information about the
universal properties of many-body systems encoded in the combinatoric
structure of the $F$ matrix. That second part requires a statistical approach
but can be efficiently computed for large system sizes by avoiding using the
$F$ matrix explicitly since its number of rows scales exponentially with the
number of levels and particles considered. It should be noted that despite the
$F$ matrix is obviously a real matrix, the renormalized energies and the
$U_{\ell}$ vectors can be complex, as we will see in the following.

In the following, we provide a framework to compute the singular decomposition
of $F$ in the case of a system of non-interacting fermions with no other
quantum number.

##  III Spinless Fermions

In the case of fermions without any additional quantum number, occupation
numbers can only be $0$ or $1$, which are the possible matrix components of
$F$. Each row of the filling matrix $F$ (i.e. some configuration of a many
body state) is a binary string. All rows share the same amount of $1$’s, to
account for the fixed number of particles $N$, there are a total of
$C_{L}^{N}$ many body state configurations and the filling matrix $F$ has
dimensions $L\times C_{L}^{N}$.

###  III.1 Right part of the SVD

We focus first on the right-singular vectors $V_{\ell}$ and the singular
values $\sigma_{\ell}$. Both can be computed by considering the
eigendecomposition of the square matrix $F^{\dagger}\cdot F$ which has
dimensions $(L\times L)$. Each matrix element $(F^{\dagger}\cdot F)_{i,j}$ is
the scalar product between the $i^{th}$ and $j^{th}$ columns of $F$, each
encoding the filling factor of the $i^{th}$ and $j^{th}$ single body states
along all of the many body states. If $i=j$, simple combinatorics
considerations indicate there are $C_{L-1}^{N-1}$ many body levels with the
$i^{th}$ single body level filled with one particle. If $i\neq j$, there are
$C_{L-2}^{N-2}$ many body levels with both the $i^{th}$ and $j^{th}$ single
body levels filled with a particle. From this reasoning we deduce that the
$F^{\dagger}\cdot F$ matrix has the form of a circulant matrix with only two
possible values :

$F^{\dagger}\cdot F=\begin{pmatrix}a&b&\dots&b\\\ b&\ddots&&\vdots\\\
\vdots&&\ddots&b\\\ b&\dots&b&a\end{pmatrix}_{L\times L}$

where $a=C_{L-1}^{N-1}$ and $b=C_{L-2}^{N-2}$. Circulant matrices are well
known and can be diagonalized using Fourier modes. Let $\omega=e^{i2\pi/L}$,
then eigenvectors of $F^{\dagger}\cdot F$, which are the right-singular
vectors $V_{\ell}$ of $F$, are given by :

$V_{\ell}=\frac{1}{\sqrt{L}}\left(1,\omega^{\ell},\omega^{2\ell},\dots,\omega^{(L-1)\ell}\right).$
(4)

The eigenvalues of $F^{\dagger}\cdot F$ take the following form :

$\lambda_{\ell}=C_{L-1}^{N-1}+C_{L-2}^{N-2}\sum_{k=1}^{L-1}\omega^{k\ell}$

from which only two distinct values can arise. After taking the square root to
obtain singular values of $F$, we get :

$\displaystyle\sigma_{0}$ $\displaystyle=\sqrt{NC_{L-1}^{N-1}}$
$\displaystyle\text{ with multiplicity }1$ $\displaystyle\sigma_{\ell\geq 1}$
$\displaystyle=\sqrt{C_{L-2}^{N-1}}$ $\displaystyle\text{ with multiplicity
}L-1$ (5)

One can easily check that the $V_{\ell}$ vectors form an orthonormal base
associated to these eigenspaces, the $\ell=0$ case matching the
$1$-dimensional eigenspace spanned by $V_{0}$ and the $\ell\geq 1$ case
matching the $L-1$-dimensional eigenspace spanned by $V_{\ell\geq 1}$.

###  III.2 Left part of the SVD

In order to compute the left-singular vectors $U_{\ell}$, we avoid dealing
with the matrix $F\cdot F^{\dagger}$ which is much larger than
$F^{\dagger}\cdot F$, and prefer using the identity $U_{\ell}=F\cdot
V_{\ell}$. We can write the components of $U_{\ell}$ as follows :

$U_{\ell}^{p}=\frac{1}{\sqrt{L}}\sum_{k=0}^{L-1}F_{k}^{p}\omega^{k\ell}$

essentially, the inner product between $V_{\ell}$ and the configuration of the
$p^{\text{th}}$ many body state.

We will now investigate two symmetries which provide a good understanding of
universal properties of $F$ and allow to reduce the very large set of
configurations to a more manageable size for numerical applications.

###  III.3 The $k$-symmetry

The first symmetry we observe is on the $V_{\ell}$ vectors containing Fourier
modes based on $L^{\text{th}}$ roots of unity. Since they contain components
of the form $\omega^{k\ell}$, when calculating the inner product of $V_{\ell}$
against a configuration represented by a binary string, a circular permutation
of those bits will multiply the Fourier modes by a power of $\omega$, i.e. it
will rotate the result $U_{\ell}^{p}$ in the complex plane by a complex factor
which is a power of $\omega^{\ell}$. We will denote the associated angle
$\theta_{\text{step}}=2\pi\ell/L$. Note that it does not depend on $N$ and
remains valid for all components of a given $U_{\ell}$. It follows that all
the components of $U_{\ell}$ lie on various circles in the complex plane, and
are located at relative angles which are multiples of $\theta_{\text{step}}$.
We call this structure the $k$-symmetry, which is illustrated on Fig. 1.

The $k$-symmetry suggests to group configurations together in equivalence
classes defined by the underlying relation of circular permutations. Each
class contains configurations that can be transformed into each other by
circular permutation, we choose the lowest of them in lexicographic order as
class representatives which we will refer to as ”seeds”. All members of such
an equivalence class, when taking the inner product with $V_{\ell}$, will
produce $U_{\ell}$ components which lie on the same circle in the complex
plane, for any fixed given $\ell$. Note that two distinct classes can still
map to the same circle by having the same radius.

There will be two kind of such equivalence classes : the non-degenerate
classes contain exactly $L$ binary strings when the bits don’t provide any
additional symmetry, while the degenerate classes contain less than $L$
strings as the bits expose a shorter pattern that repeats under circular
permutation, for instance $(0,0,1,1,0,0,1,1)$ is an $8$-bits sequence but
loops back early with a circular permutation by $4$ positions. This
distinction will become important later on.

Next, we notice that equivalence classes are closely related to the factors of
$L$. In particular, if $L$ is a prime number, all classes are non-degenerate.
Degenerate classes can only have a cardinal which has a common factor with
$L$. This can be deduced from the action of circular permutations on strings
of various sizes. This is of particular physical meaning : the prime
decomposition of $L$ is a decisive feature for the symmetries of the filling
matrix and its singular value decomposition.

Finally, a counting argument can be made about those equivalence classes by
using the Pólya enumeration theorem, from which we can estimate the number of
equivalence classes to be of the order of $C_{L}^{N}/L$.

###  III.4 The $\ell$-symmetry

On the other hand, we can look at what happens when we consider a different
Fourier mode and go from some $\ell$ to another $\ell^{\prime}$. As $\omega$
is by definition a primitive $L^{th}$ root of unity (i.e. $\omega^{L}=1$ but
$\omega^{p}\neq 1$ for all $p<L$), then each $L^{th}$ root of unity is a
distinct power of $\omega$. It is well known that $\omega^{\ell}$ is a
$r^{th}$ primitive root of unity for $r=L/\text{gcd}(L,\ell)$. Considering
another value $\ell^{\prime}$, if

$\gcd(L,\ell)=\gcd(L,\ell^{\prime})$ (6)

then $\omega^{\ell^{\prime}}$ is also a $r^{th}$ primitive root of unity. In
other words, each power of $\omega^{\ell}$ is in a one to one correspondance
with a power of $\omega^{\ell^{\prime}}$, so the components of
$V_{\ell^{\prime}}$ are just a permutation of the components of $V_{\ell}$. In
this case, since the permutation of $V_{\ell}$ components is equivalent to the
permutation of configurations in the $F$ matrix, then vectors $U_{\ell}$ and
$U_{\ell^{\prime}}$ are also the same up to a permutation, i.e. the
distributions of their components are identical and the resulting spectra are
identical. We call this property the $\ell$-symmetry, verified when (6) is
valid, which is illustrated in Fig.6 where all principal spectra are provided
for $L=20$ and $N=10$.

In order to proceed further, we need to introduce the language of
compositions. A composition of an integer is similar to a partition, but with
order taken into account Knuth (2015). Given some integer $n$ (in our case the
number of particles), we decompose $n$ into a sum of $p$ non-zero integer
called the ”parts”. Any such sequence is called a $p$-composition of $n$. If
we allow some of the $p$ parts to be zero, then it is called a weak
$p$-composition of $n$. Finally, if we also impose some integer $k$ such that
any part can only have a maximum value of $k$, then such list of integers is
called a $k$-restricted weak $p$-composition of $n$.

The $\ell$-symmetry simplifies the study of the principal spectra, as only a
small amount of $\ell$ values are required to obtain the complete set of
$U_{\ell}$, namely the divisors of $L$. To implement this, we can change our
point of view on binary strings, and consider them as $1$-restricted weak
$L$-compositions of the integer $N$. We then proceed by considering only the
$\ell$’s which divide $L$, and introduce the integer $q=L/\ell$. We subdivide
each binary string into $\ell$ sections of size $q$, and add them together
component by component, effectively folding them into a vector of size $q$. We
obtain $\ell$-restricted weak $q$-compositions of the integer $N$. This leads
to the following simplification :

$U_{\ell}^{p}=\frac{1}{\sqrt{L}}\sum_{k=0}^{q-1}\left(\sum_{s=0}^{\ell-1}F_{sq+k}^{p}\right)e^{i2\pi\frac{k}{q}}$

from which we redefine effective vectors $V^{\prime}_{\ell}$ of size $q$,
which now contain Fourier modes based on the $q^{\text{th}}$ roots of unity
$\omega^{\prime}=e^{i2\pi/q}$. Its components are of the form $\omega^{\prime
k}$ with $k\in\\{0,q-1\\}$.

Once again, we see that the divisors of $L$, i.e. its prime decomposition,
plays a central role. In particular, if $L$ is prime, we only need to compute
the $\ell=1$ case where the $\ell$-symmetry is trivial and the
$\ell$-restricted weak $q$-compositions are the binary strings themselves. For
a composite $L$, there is an interesting feature which can be observed from
Fig.2, where increasing $\ell$ values lead to a fast decreasing density of
$U_{\ell}$ components in the complex plane, while the occurence counts of each
points necessarily goes up : we transition from a spread distribution to a
clustered distribution, i.e. with larger degeneracies, and a lattice pattern
can emerge. The full understanding of the rich and qualitatively different
behaviors displayed on Fig.2 as $l$ increases requires further study and will
be the subject of another article.

Fig. 1:  First principal spectrum $U_{1}$ of the filling matrix $F$ for $L=8$
single body levels and $N=4$ spinless fermions. Each point is a distinct value
of the distribution of the components of $U_{1}$ in the complex plane. Circles
relate components belonging to one or more equivalence classes having the same
modulus. The number of occurences of each values are : $6$ (dark purple) for
all degenerate classes with radius $0$, $1$ (yellow) which have $1$ non-
degenerate class per circle, and $2$ (dark orange) which have $2$ non-
degenerate classes per circle.  (a)(b) Fig. 2:  Principal spectra of the
filling matrix $F$ in the case $L=20$ at half-filling $N=10$ (a) or quarter-
filling $N=5$ (b). Each graph represents components of the $U_{\ell}$ vector
in the complex plane and their occurrence numbers. All relevant $\ell$ values
with distinct greatest common divisors with $L$ are shown, other $\ell$ values
yield identical distributions through $\ell$-symmetry, according to the
criteria from Eq.(6). The full set of $\ell$ values can be found in Fig. 6.
The $\ell=0$ case is also shown and corresponds to the $1$-dimensional
eigenspace associated to the singular decomposition of $F$. One can notice
that increasing values of $\ell$ yield more clustered distributions.

###  III.5 Enumeration and statistics

Both symmetries described above simplify the study of the problem. Going from
the set of $C_{L}^{N}$ binary strings to a set of $\ell$-restricted weak
$q$-compositions is not an injection, i.e. distinct binary strings can have
the same weak composition, so occurrence numbers should be tracked. Let’s
denote by $\\{m_{0},\dots,m_{q-1}\\}$ the parts of an $\ell$-restricted weak
$q$-composition, with $m_{k}\in\\{0,\ell\\}$ and their sum adding up to $N$.
Then, the total number of binary strings mapped to the same composition is
given by a winding factor :

$Q=\prod_{k=0}^{q-1}C_{\ell}^{m_{k}}$ (7)

From this point on, all ”folded” configurations are now represented by a
choice for each of the $m_{k}$ values and we use $V^{\prime}_{\ell}$ vectors
of size $q$ containing $q$-roots of unity $\omega^{\prime k}$ for
$k\in\\{0,q-1\\}$. For notation simplicity, we will now drop the primes and
simply redefine $V_{\ell}$ and $\omega$.

The next step is to apply the $k$-symmetry on top of compositions, by again
noticing that any circular permutation of a composition is only a
multiplication of $V_{\ell}$ by a power of $\omega$ : we group compositions
together in equivalence classes related by circular permutations. This is
where an interesting feature comes in : by applying the $k$-symmetry on
compositions, it can be shown that configurations inside degenerate classes
will always have a null projection on the Fourier modes of $V_{\ell}$. We only
need to count them to obtain the occurence count of zero components in
$U_{\ell}$.

We now have all the tools required for our framework. We can describe a simple
recipe to obtain the exact distribution of all $U_{\ell}$ vectors. First, we
look at the value of $L$ and construct the list of all its divisors : this
gives us the $\ell$ values we have to deal with. For each of them, we compute
$q=L/\ell$ and proceed by building the list of seeds : all $\ell$-restricted
weak $q$-compositions which are not related by circular permutation. For each
seed, we first look at the size of its equivalence class. If it is degenerate,
i.e. lower than $q$, then the class lies on a $0$-radius ”circle” which is a
point at origin in the complex plane. In this case, we only need to recover
the occurrence count which is the winding factor $Q$ from Eq.(7) times the
degenerate class size. Otherwise, the class is non-degenerate and has size $q$
: the total occurrence count is $Q$ times $q$, each vertex having an
occurrence count of $Q$. We also compute the inner product of the seed with
$V_{\ell}$ directly and extract the module and argument of the result. The
module needs to be normalized by $1/\sqrt{\ell}$ to give us the radius of this
seed’s circle. The argument gives us one point on the circle from which we can
reconstruct all other points coming from this seed’s equivalence class by
using increments of $\theta_{\text{step}}$, i.e. the $k$-symmetry. We now know
everything about that particular seed and its class. While processing all
unique seeds and classes, we can accumulate the ones found to lie on the same
circle, i.e. having the same radius, and simply add up occurrence counts
accordingly. The end result is a list of all the circles, their respective
radius, angular positions of all vertices if needed, and occurrence numbers
(by point or by circle, whichever is needed). This gives the exact
distribution of $U_{\ell}$ components. Finally, all other $\ell$ values which
were skipped are recovered from $\ell$-symmetry by the criteria from Eq.(6).

As an example, let us consider the case $L=8$ at half-filling $N=4$. There are
a total of $70$ possible many body configurations written as binary strings.
The divisors of $L$ are $\\{1,2,4\\}$. For $\ell=1$ (i.e. binary strings), the
$k$-symmetry gives $8$ non-degenerate classes of size $8$ and $2$ degenerate
classes of size $4$ and $2$, for a total of $10$ seeds. For $\ell=2$, applying
both symmetries, we have $2$-restricted weak $4$-compositions, $4$ non-
degenerate classes of size $4$ and $2$ degenerate classes of size $2$ and $1$,
for a total of $6$ seeds. For $\ell=4$, applying both symmetries, we have
$4$-restricted weak $2$-compositions, $2$ non-degenerate classes of size $2$
and $1$ degenerate class of size $1$, for a total of $3$ seeds. As one can
see, the number of configurations to consider is significantly reduced. For
larger systems, e.g. $L=20$ at half-filling $N=10$, the number of many body
configurations written as binary strings is $184756$. In the worse case,
$\ell=1$, this is reduced to $9252$ binary seeds. In better cases like
$\ell=4$, using $\ell$-symmetry leads to a much smaller set of $381$
compositions which is further reduced to only $77$ seeds.

For even larger systems, ie. $L$ of the order of $100$, the lower $\ell$
values still require significant computational time, mostly for $\ell=1$,
unless they are truncated or a statistical approach is added to sample the set
of seeds. However, it should be noted that those computations are independent
of the single body spectrum: they should be performed only once for relevant
values of $L$ and $N$ (ie. half-filling, quarter-filling, …) and re-used for
many different systems. Once the distribution of $U_{\ell}$ components is
known, the many-body spectrum and its DoS are obtained from Eq. (3). The
spectrum for $l=1$ which is of particular importance for the applications is
displayed on Fig.5.

##  IV Applications

A variety of systems can be described in terms of spinless fermions, including
hard core bosons Girardeau (1960), Mott insulators in ladders Donohue  _et
al._ (2001), spin liquids Imambekov  _et al._ (2012), and, very recently, they
served as the basis to study topological phases Turner  _et al._ (2011) and
systems supporting Majorana fermions Kitaev (2001); Alicea (2010).

To illustrate our method, we will apply it to two classes of problems: (i)
tight-binding describing electrons in a one-dimensional chain and in a square
lattice, and (ii) the transverse Ising field chain that describes a critical
1D spin chain, and which can be mapped into a single-particle problem using
the Jordan Wigner (JW) transformation. In both, we consider periodic boundary
conditions (PBC).

In the case (i), the Hamiltonian reads

$\displaystyle H=-t\sum_{\braket{i,j}}(c_{j}^{\dagger}c_{j}+\text{H.c.}),$ (8)

where the sum $\braket{i,j}$ runs over first-neighbor sites,
$c_{i}/c_{i}^{\dagger}$ are the annihilation/creation operators, and $t$ is
the hopping amplitude. In 1D, the single-body energies are

$\displaystyle\epsilon_{k}$ $\displaystyle=-2t\cos(k),$ (9)

the momenta $k$ depends on the boundary conditions of the model. For periodic
boundary conditions (PBC) $k=\frac{2\pi n}{L}$, with $n=0,...,L-1$.

The trigonometric form of the dispersion relation associated with the fact
that the matrix Hamiltonian is circulant has an interesting implication to the
renormalized energies in Eq.(3). The periodicity shared between the single-
body energies and the Fourier modes allows for the cancellation of all
renormalized energies $\tilde{\epsilon}_{\ell}$ except those associated with
$\ell=1$ and $\ell=L-1$, which are symmetric to each other. We can show that
these surviving contributions are equal to
$\tilde{\epsilon}_{1}=\tilde{\epsilon}_{L-1}=-t\sqrt{L}$. As a consequence,
the Many Body DoS involves only one principal spectrum, the $l=1$ and is
displayed on Fig. 5.

In 2D, the band structure contains $L^{2}$ single-body energies given by

$\displaystyle\epsilon_{k_{x},k_{y}}$
$\displaystyle=-2t\cos(k_{x})-2t\cos(k_{y}),$ (10)

where the momenta in each direction are $k_{x}=\frac{2\pi n_{x}}{L}$ and
$k_{y}=\frac{2\pi n_{y}}{L}$, with $n_{x},n_{y}=0,...,L-1$.

Similarly to the one-dimensional case, the renormalized band structure of the
square lattice will be non zero for only a few $\ell$’s, due to the high
degeneracy of the single-body energies. In particular,
$\epsilon_{k_{x},k_{y}}=0$ has degenegeracy $2(L-1)$, whereas the only non-
degenegerate energies occur at the bottom and at the top of the band, where
$\epsilon_{k_{x},k_{y}}=-4J$ and $\epsilon_{k_{x},k_{y}}=4J$, respectively. By
flattening the 2D band structure as a 1D vector with $L^{2}$ entries, one can
show that the $\tilde{\epsilon}_{\ell}$ do not vanish for
$\ell=nL+1,n=0,1,...,L/2-1$ and $\ell=mL-1,m=1,2,...,L/2$, and for $\ell=L$,
which is a special case in which $\tilde{\epsilon}_{\ell}$ is real. The band
structure in Eq. (10) and the non-vanishing renormalized energies
$\tilde{\epsilon}_{\ell}$ for a square lattice with $L=100$ sites is shown in
Fig. 3. The real part decays fast with $\ell$, while the imaginary part
converges to a constant, with oscillations decreasing with $\ell$.

Fig. 3:  Band structure and renormalized energies of a homogeneous square
lattice with $L=100$ sites in each direction. Top panels show the single-body
band structure $\epsilon_{k_{x},k_{y}}$ as a surface in 3D (left) and as a
contourplot (right). Bottom panels display the real and imaginary part of the
renormalized energies $\tilde{\epsilon}_{\ell}$. Note that only $2(L+1)$
renormalized energies indexed as $\ell=nL\pm 1$ do not vanish.

In the case (ii), the Hamiltonian is

$\displaystyle H$
$\displaystyle=-J\sum_{j}\sigma_{j}^{x}\sigma_{j+1}^{x}-h\sum_{j}\sigma_{j}^{z},$
(11)

where $\sigma^{d}$ $d=x,y,z$ are the Pauli matrices, $J$ is the coupling and
$h$ is the transverse field.

The critical point occurs at $h=J$. After employing the JW transformation, we
can obtain the dispersion relation of this model as follows

$\displaystyle\epsilon_{k\pm}$ $\displaystyle=\pm
2J\sqrt{h^{2}+1-2h\cos(k)}-2Jh,$ (12)

with $k=\frac{2n\pi}{L}$ and $n=0,...,L-1$ being the momentum. Note that now,
two bands, one positive and the other negative, contribute to the density of
states.

In this case, the structure of the momenta $k$ corresponding to the Fourier
modes result in most of the renormalized energies to be non-zero, however only
odd indexes contribute. For the real part, all
$\text{Re}[\tilde{\epsilon}_{\ell}]=4/\sqrt{2L}|h/J-1|$ with odd $\ell$ are
independent of $\ell$. The factor in modulus introduces a symmetry around the
critical point $h_{c}=J$, so that
$\tilde{\epsilon_{\ell}}[h-h_{c}<0]=\tilde{\epsilon_{\ell}}[h-h_{c}>0]$. For
the imaginary part, the $\text{Im}[\tilde{\epsilon}_{\ell}]$ are non-zero but
decay very fast. The decay rate is amplified as the system size increases. At
the critical point $h=J$, the only surviving single-body energies are
$\tilde{\epsilon}_{1}=-\tilde{\epsilon}_{L-1}=2i\sqrt{2L}$. In other words, we
recover the same renormalized single body energies as for the 1D tight binding
model.

This can be observed in Fig. 4, panel (b), where we show the real and
imaginary parts of $\tilde{\epsilon}_{\ell}$ for $h/J\in[0.5,1.5]$.

$\epsilon_{k}$ Fig. 4:  Band structure and renormalized single-particle
energies $\tilde{\epsilon}_{\ell}$ for an Ising chain under a transverse field
with strength $h$. The number of sites/levels is $L=100$. Top contour shows
the band dispersion as a function of $k$ and $h$. Bottom panels display the
real and imaginary parts of the $\tilde{\epsilon}_{\ell}$’s for $\ell$ odd.
All components with $\ell$ even vanish. The imaginary part decays as a
function of $\ell$. Both real and imaginary parts are degenerate around
$h_{c}=J$, i.e.
$\tilde{\epsilon_{\ell}}[h-h_{c}<0]=\tilde{\epsilon_{\ell}}[h-h_{c}>0]$. At
the critical point, the dominant renormalized energy is
$\tilde{\epsilon}_{1}$, all other renormalized energies can be neglected in a
first approximation. As a result, the MBDoS will be similar to the one of the
1D tight binding model.

Using these previous results, we have all ingredients to compute the MBDoS for
the applications aforementioned. In particular, for the 1D tight binding model
and the transverse Ising chain at the critical point which have the same
renormalized single body energies. As a result, they have the same MBDoS
displayed on Fig. 5.

(a)(b) Fig. 5:  Many-body DoS for a 1D tight-binding chain with $L=20$ sites
at half-filling $N=10$ (a), or quarter-filling $N=5$ (b). These histograms of
the many body energies are are obtained by applying Eq.(3) where the only non-
zero renormalized energies $\tilde{\epsilon}_{\ell}$ are for $\ell=1$ and
$\ell=L-1$ (see Sec. IV) and with $80$ bins from the lowest to the highest
many body energies. The MBDoS is then normalized to a probability density.
Note that the particle-hole symmetry explains the even parity of the profile.

##  V Conclusion

In the present paper, we explored a novel approach to compute the many-body
density of states of quantum systems whose Hamiltonian can be mapped into non-
interacting spinless fermions.

We show that the many-body spectrum can be expanded over the principal
components of a filling matrix encoding the allowed many body states. These
principal components describe spectral properties of systems of spinless
fermions only depending on the numbers of particles and one body levels. This
new spectral decomposition of many body states is weighted by a renormalized
single body band structure, which acts as a filter for relevant energy scales.

For gapless systems, such as the tight-binding and the critical transverse
field Ising chains, we demonstrated that only two renormalized energies are
non zero. Even in more general scenari, such as the square lattice or the
Ising chain away from the critical point, many renormalized energies still
vanish. In all cases, this will significantly reduces the number of relevant
spectral components of the filling matrix involved in the calculation of the
many body density of states.

Our framework can be extended to include additional quantum numbers like spin,
and to handle bosonic systems.

##  VI Acknowledgements

We acknowledge support from the Leverhulme Trust under grant RPG-2020-094.

(a)(b) Fig. 6:  Full spectrum of components of $U_{\ell}$’s for $L=20$ at
half-filling $N=10$. All $\ell$ values with the same greatest common divisor
with $L$ are associated to the same principal spectrum $V_{\ell}$,
illustrating the $\ell$-symmetry discussed in Sec.III.

## References

  * Zelevinsky and Horoi (2019) V. Zelevinsky and M. Horoi, Progress in Particle and Nuclear Physics 105, 180 (2019).
  * Ashcroft  _et al._ (1976) N. Ashcroft, A. W, N. Mermin, N. Mermin, and B. P. Company,  _Solid State Physics_, HRW international editions (Holt, Rinehart and Winston, 1976).
  * Gogolin and Eisert (2016) C. Gogolin and J. Eisert, Reports on Progress in Physics 79, 056001 (2016).
  * Santos and Rigol (2010) L. F. Santos and M. Rigol, Phys. Rev. E 81, 036206 (2010).
  * Borgonovi  _et al._ (2016) F. Borgonovi, F. Izrailev, L. Santos, and V. Zelevinsky, Physics Reports 626, 1 (2016), quantum chaos and thermalization in isolated systems of interacting particles.
  * Polkovnikov  _et al._ (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, Annual Review of Condensed Matter Physics 6, 15 (2015), https://doi.org/10.1146/annurev-conmatphys-031214-014726 .
  * Abanin  _et al._ (2019) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * Ithier and Benaych-Georges (2017) G. Ithier and F. Benaych-Georges, Phys. Rev. A 96, 012108 (2017).
  * Ithier  _et al._ (2017) G. Ithier, S. Ascroft, and F. Benaych-Georges, Phys. Rev. E 96, 060102 (2017).
  * Weiße and Fehske (2008) A. Weiße and H. Fehske, in  _Computational many-particle physics_ (Springer, 2008) pp. 529–544.
  * Weiße  _et al._ (2006) A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, Rev. Mod. Phys. 78, 275 (2006).
  * Bethe (1936) H. A. Bethe, Phys. Rev. 50, 332 (1936).
  * Volya and Zelevinsky (2006) A. Volya and V. Zelevinsky, Phys. Rev. C 74, 064314 (2006).
  * Hillman and Grover (1969) M. Hillman and J. R. Grover, Physical Review 185, 1303 (1969).
  * Berger and Martinot (1974) J. Berger and M. Martinot, Nuclear Physics A 226, 391 (1974).
  * Jacquemin and Kataria (1986) C. Jacquemin and S. Kataria, Zeitschrift für Physik A Atomic Nuclei 324, 261 (1986).
  * Bohr and Mottelson (1998) A. N. Bohr and B. R. Mottelson,  _Nuclear Structure (in 2 volumes)_ (World Scientific Publishing Company, 1998).
  * Knuth (2015) D. E. Knuth,  _The art of computer programming, Volume 4_ (Addison-Wesley Professional, 2015).
  * Girardeau (1960) M. Girardeau, Journal of Mathematical Physics 1, 516 (1960), https://doi.org/10.1063/1.1703687 .
  * Donohue  _et al._ (2001) P. Donohue, M. Tsuchiizu, T. Giamarchi, and Y. Suzumura, Phys. Rev. B 63, 045121 (2001).
  * Imambekov  _et al._ (2012) A. Imambekov, T. L. Schmidt, and L. I. Glazman, Rev. Mod. Phys. 84, 1253 (2012).
  * Turner  _et al._ (2011) A. M. Turner, F. Pollmann, and E. Berg, Phys. Rev. B 83, 075102 (2011).
  * Kitaev (2001) A. Y. Kitaev, Physics-Uspekhi 44, 131 (2001).
  * Alicea (2010) J. Alicea, Phys. Rev. B 81, 125318 (2010).