ampmtime

# Consequences of Multiple Axions in Theories with Dark Yang-Mills Groups

Manuel Ettengruber<EMAIL_ADDRESS>Arnold Sommerfeld Center, Ludwig-
Maximilians-Universität, Theresienstraße 37, 80333 München, Germany,  Max-
Planck-Institut für Physik, Föhringer Ring 6, 80805 München, Germany
Emmanouil Koutsangelas<EMAIL_ADDRESS>Arnold Sommerfeld Center, Ludwig-
Maximilians-Universität, Theresienstraße 37, 80333 München, Germany,  Max-
Planck-Institut für Physik, Föhringer Ring 6, 80805 München, Germany

(, )

###### Abstract

General consistency requirements of Quantum Gravity demand the existence of
one axion per Yang-Mills group. In this work, we consider theories with dark
Yang-Mills sectors and investigate general phenomenological implications of
these necessary axions. We carry out computations for two simple models,
namely a pure Yang-Mills sector and $N$ exact Standard Model copies. For the
former, the misalignment mechanism results in a minimal dark confinement scale
$\Lambda_{\rm conf}\gtrsim 1\,\mathrm{eV}$ if the dark sector axion is
supposed to make up the dark matter. For the latter, the misalignment
mechanism without fine-tuning of the initial misalignment angle places an
upper bound on $N$ below the species bound. When the PQ symmetries are broken
during inflation, the collective isocurvature fluctuations do not necessarily
tighten the bound on the inflationary Hubble scale arising from a single
axion. We also point out that axion stars collectively made from axions of
different dark sectors with a suppressed mass spectrum are not possible.
Lastly, for the two models at hand, intersector interaction through axion
kinetic mixing leads to the existence of two distinct axion states. For a
single dark YM sector, the upper bound $\Lambda_{\rm conf}\lesssim
10^{12}\,\mathrm{GeV}$ emerges from the stability requirement of the dark
sector axion. For $N$ exact SM copies, the mass and photon coupling of the
second state is completely determined after a potential measurement of the
analogous parameters of the first axion.

###### Contents

  1. I Introduction
  2. II Cosmological Implications
    1. II.1 Cosmological Framework
    2. II.2 Bounds on the number of species
    3. II.3 The Axion Mechanism
    4. II.4 Misalignment: Exact SM Copies
    5. II.5 Misalignment: Pure YM Dark Sector
    6. II.6 Isocurvature Perturbations
    7. II.7 $N$-Boson Stars
  3. III Kinetic Mixing between Axions
    1. III.1 Kinetic Mixing as Intersector Interaction
    2. III.2 Modification of Axion Physics for SM Copies
    3. III.3 Modification of Axion Physics for one YM Sector
  4. IV Conclusions

##  I Introduction

Attempted measurements of the neutron electric dipole moment indicate no CP
violation in the strong interaction. Since QCD includes CP violating processes
proportional to a parameter $\theta$, the non-observation results in the bound
$\theta<10^{-10}$ ThetaBound . This inexplicable smallness is referred to as
the strong CP problem. While it looks like a naturalness problem this is
actually misleading. It would be a naturalness problem if fine-tuning was
necessary due to unacceptably large quantum corrections (like for the Higgs
mass). For $\theta$ the quantum corrections from the Standard Model are much
smaller than the measured bound ThetaCorrections , so that from the point of
view of naturalness it is rather a small value puzzle than a problem.

However, the situation dramatically changes when quantum gravity is taken into
account dvali2022strong . In recent years, evidence has been provided that
quantum gravity is incompatible with (meta-)stable de Sitter vacua
Dvali:2013eja ; Dvali:2018jhn . In particular, this is organic in the S-matrix
formulation of quantum gravity, which is only compatible with asymptotically
flat space-times Dvali:2020etd . This has fundamental implications for axion
physics Dvali:2018dce ; dvali2022strong . The point is that in its essence
$\theta$ is not just a parameter but a label for QCD vacua. These vacua are
non-degenerate. In particular, according to the Vafa-Witten theorem
PhysRevLett.53.535  the global minimum is the CP conserving one at $\theta=0$.
Other vacua have higher energies. Correspondingly, as argued in Dvali:2018dce
; dvali2022strong , if all theta vacua were to be physical, their non-
degeneracy would imply the existence of de Sitter vacua among them. Moreover
dvali2022strong , even at the expense of arbitrary fine-tuning, only one out
of infinitely many vacua could be made viable for the S-matrix formulation.
Others would end up in asymptotically non-flat cosmologies. Correspondingly,
it was concluded in Dvali:2018dce ; dvali2022strong  that $\theta$ vacua must
be unphysical in gravity.

This means that when gravity is included the strong CP problem becomes a
consistency problem and the non-degenerate vacuum structure must be removed,
leaving only the viable vacuum. Therefore, a mechanism to select the right
vacuum becomes indispensable Dvali:2018dce ; dvali2022strong , which was
provided by Peccei and Quinn (PQ) PQMechanism ; PQMechanism2 . By introducing
a non-linearly realized $\mathrm{U}(1)_{\rm PQ}$ that is anomalous with
respect to QCD, the $\theta$-parameter is dynamically set to zero by the
pseudo-Goldstone boson of $\mathrm{U}(1)_{\rm PQ}$ also known as the axion
WeinbergAxion ; Wilczek:1977pj .

Viewed in this light, the strong CP problem is not exclusive to QCD but
extends to any non-Abelian Yang-Mills (YM) group. Every YM group includes a
vacuum angle that leads to an unacceptable de Sitter-type vacuum structure and
must therefore be removed dvali2022strong .

Using the PQ mechanism to achieve this, requires one axion per YM group
because a single axion cannot set multiple $\theta$-parameters to zero. In
particular, this is the case even when the different YM sectors are related by
a discrete symmetry dvali2005threeform , e.g. an exact permutation symmetry
among $N$ Standard Model copies Dvali:2007iv . This is due to the fact that
the topological susceptibilities in different YM sectors break the discrete
symmetry spontaneously. Correspondingly, a single axion can only eliminate one
combination of the vacuum angles. The other combinations will remain physical
and unconstrained, resulting in two problems. First, even ignoring gravity,
this would imply that a single axion cannot solve the small $\theta$ puzzle
for multiple YM groups dvali2005threeform ; Dvali:2007iv . Secondly, allowing
these vacuum angles to persist, would necessarily introduce the asymptotically
non-flat vacua that would conflict with the requirement of an S-matrix
formulation, hence requiring one axion per YM group dvali2022strong .

In this work, we consider theories in which dark matter includes one or more
YM sectors and investigate the generic phenomenological consequences of the
necessary axions. To our knowledge, axions have not been taken into account
due to gravitational consistency reasons, which impose certain properties on
them.

To be precise, we consider models with $N$ sectors $\mathcal{L}_{i}$, each
including at least one YM subgroup. Without loss of generality, we choose our
sector to be labeled by $i=1$. All the others collectively make up the dark
matter and we will refer to them as hidden or dark sectors. The sectors
communicate via gravity and additional non-gravitational interactions, which
we encode in $\mathcal{L}_{\rm mix}$. The associated Lagrangian can be written
as

$\mathcal{L}=\sum_{i=1}^{N}\mathcal{L}_{i}+\mathcal{L}_{\rm mix}\;.$ (1)

The possible models can be divided in two categories. The first is a framework
in which the dark sector consists of $N$ copies of the standard model related
by an exact permutation symmetry Dvali:2009fw ; Dvali:2009ne . This case is
most restrictive due to exact symmetry. Large values of $N$ are motivated by
the Hierarchy problem Dvali:2007hz ; Dvali:2007wp  but we shall also cover the
case with smaller $N$. Such were discussed in different contexts, e.g., in
Blinnikov:1982eh ; Kolb:1985bf .

The second category includes all other possibilities. The dark sector in such
cases is unrelated to the SM and can contain very different gauge structures.
Such cases are expected to be rather generic in string theory
compactifications PhysRevLett.54.502 ; Dixon:1985jw ; Lebedev:2006kn ;
Braun:2005ux .

Before proceeding let us note that there exists a universal constraint on the
extent of dark sectors in the form of the bound on the total number of
particle species, $N_{s}<10^{32}$ Dvali:2007hz ; Dvali:2007wp . This bound is
fully non-perturbative in its nature and originates from black hole physics,
which implies that in theory with $N_{s}$ particle species the quantum gravity
scale is lowered to $M_{*}=M_{P}/\sqrt{N_{s}}$ Dvali:2007hz . From current
phenomenology, this translates as the above bound on $N_{s}$, since in this
extreme case the scale would be lowered to TeV. This would solve the hierarchy
problem, which provides the main motivation for dark sectors with many
species. Another important non-perturbative bound on the number of species
originates from the constraints on the species entropy Dvali:2020wqi . A more
detailed discussion on both bounds on the number of species will be given
later.

For concreteness, we will look at two simple models: One pure YM sector and
$N$ exact SM copies. The former has been used to model dark matter due to the
dynamics of YM theories, which for instance may result in a potential phase
transition that produces stochastic gravitational waves during the early
universe Halverson:2020xpg ; PhysRevD.104.035005 . The latter was originally
motivated by the many-species solution to the hierarchy problem Dvali:2007hz ;
Dvali:2007wp . In dealing with extended dark sectors, one can distinguish
between universal effect and model-dependent features characteristic of
particular scenarios. In particular, a universal feature is the modification
of black hole physics at scales below the species scale $M_{*}$ (see,
Dvali:2008fd ). In addition, there are model dependent consequences. When
applied to model the dark matter Dvali:2009fw , it gives rise to many
interesting phenomena such as neutron and neutrino oscillations Dvali:2009ne ;
Ettengruber:2022pxf , modified black hole physics
MicroBlackHolesDemocraticTransition , and compact dark matter objects
CompactDarkMatterObjectsViaN . Even though we perform calculations for these
two models, we shall try to keep our discussion maximally general and to point
out the model independent predictions.

By adding to the dark YM sectors the necessary axions, we find several
phenomenological consequences such as constraints on the parameter space of
the dark sectors, potential new small structures, and new experimental
signatures that we briefly summarize here.

The first consequence arises due to axion production via the misalignment
mechanism MisalignmentPRESKILL ; MisalignmentDINE ; MisalignmentAbbott .
Today’s axion energy density can of course be less than the observed dark
matter density but it should not be higher. An important parameter is the
initial misalignment angle of the axion field at the onset of the latest
oscillations. As shown in Dvali:1995ce  this quantity is highly sensitive to
earlier history and can be made arbitrarily small due to the early relaxation
mechanism. We shall therefore treat it as a free input parameter for each
axion sector. In the absence of an early relaxation mechanism Dvali:1995ce ,
we find the following. For a single dark YM sector, this requirement puts a
lower bound to the confinement scale of the dark YM sector, i.e. $\Lambda_{\rm
conf}\gtrsim\,\mathrm{eV}$, if the dark sector axion is required to be the
dominant contribution to the dark matter. In the case of more YM sectors the
production potentially happens in several sectors. For $N$ exact SM copies,
since each SM includes a QCD subgroup, this results in an upper bound on $N$.
Notably, the misalignment mechanism is independent of the thermalization of
the dark sectors and thus applies even for large values of $N$ when the
sectors are very dilute.

The presence of additional axions has another cosmological impact when the
axions are present during inflation, massless during this time, and not
driving inflation. Any axion that fulfills these conditions will be subject to
quantum fluctuations of isocurvature-type. Since isocurvature fluctuations
leave a special imprint in the CMB, their collective influence could tighten
the bound on the inflationary Hubble scale arising from a single axion. We
show that this is not necessarily the case, meaning that we identify regions
in the parameter space where the isocurvature bound is avoided.

In models with several dark sectors, it is possible that particles from
different sectors collectively form structures CompactDarkMatterObjectsViaN .
While in CompactDarkMatterObjectsViaN  this phenomenon was studied with SM-
like dark sectors, it does also apply to other compositions of the dark
sectors. We show that this phenomenon can result in boson stars that are
collectively made from scalars of different dark sectors. In accordance with
the literature, we will call these $N$-boson stars. The interesting feature of
such objects is that their size and mass are smaller by a factor of
$1/\sqrt{N}$ compared to the single sector case. However, for axions, this
phenomenon does not occur due to attractive self-interactions.

Lastly, we ask the question if non-gravitational communication between the
dark sectors can induce a modification of the axion physics in our sector. A
potential source for this is kinetic mixing between axions of different
sectors. We study this in the context of $N$ exact SM copies where the
underlying discrete symmetry forces all kinetic mixing parameters to be equal.
We argue that kinetic mixing emerges in both, the standard PQ scenario as well
as in the two-form gauge implementation of the axion, and that in the standard
PQ scenario, the parameters encoding the strength of the kinetic mixing,
$\epsilon$, are suppressed by $1/N$ due to perturbative unitarity. The kinetic
mixing results in $N-1$ degenerate states and one special state, which is
lighter by a factor of $1/\sqrt{2}$ and has couplings suppressed by
$1/\sqrt{N}$ (calculated for large values of $N$) as compared to the heavier
axion states. Interestingly, the lighter state is completely determined by the
heavier states, meaning that its mass and couplings are predicted once the
mass and couplings of the heavier states are known. For a single dark YM
sector, the mixing induces a decay channel of the dark sector axion into
photons. If the dark sector axion is supposed to make up the dark matter, its
stability requirement results in an upper bound on the dark confinement scale
of $\Lambda_{\rm conf}\lesssim 10^{12}\,\mathrm{GeV}$.

It is also worth mentioning that in the case of many dark sectors including
their associated axions, we arrive at a situation that is reminiscent of the
string theory inspired axiverse Arvanitaki:2009fg . Interestingly, the axions
in our framework arise from a motivation that is completely different from the
one in string theory. The phenomenology of several axions in the context of
the axiverse was for instance studied in Mack:2009hs .

The paper is organized as follows. To begin with, in Sec. II we start by
reviewing the cosmological framework of $N$ dark sectors. We then discuss how
to implement the PQ mechanism in each sector and show the bounds arising from
the misalignment production and isocurvature perturbations of $N$ axions.
Lastly, we discuss collective bound states of multiple axions. In Sec. III.1
we study the consequences of kinetic mixing between $N$ axions and point out
the differences for different types of experiments. Finally, in Sec. IV we
summarize our results and give an outlook.

##  II Cosmological Implications

###  II.1 Cosmological Framework

To begin with, let us clarify the behaviour of dark sectors within a
cosmological framework. Let us, for now, neglect intersector interactions
except gravity by setting $\mathcal{L}_{\rm mix}=0$ and discuss them
separately in Sec. III.1. This will not affect the conclusions of this
section.

For the viability of these models in a cosmological framework, the dark
sectors must behave as a cold and pressureless fluid that forms stable dark
matter halos. In addition, if the dark sectors include any massless fields,
they need to avoid the bound on the number of massless particles during
nucleosynthesis. These conditions might seem difficult to realize, especially
in the presence of non-trivial interactions in the dark sectors. We will apply
the solution proposed in Dvali:2009fw  to achieve this.

There, inflation was used to populate and reheat the sectors differently. The
crucial point is that the reheating into the dark sectors is far less
efficient than in the visible sector. The consequent low reheating temperature
in the dark sectors is used to avoid the nucleosynthesis bound, suppress
cooling processes that result in the collapse of dark matter halos, and
generally result in a different behavior throughout the cosmological history
despite the potential similarity to the visible sector. For large $N$, when
each sector is so dilute that it is problematic to regard them as thermalized
at all, the behaviour as a pressureless fluid arises from kinematics, i.e.
from the fact that particles from the same sector barely meet each other.

The situation is unchanged even when starting with equal sectors. Cosmological
evolution breaks the symmetry in a sense by resulting in different energy
densities and temperatures. The field theoretic parameters of each sector,
however, remain equal.

###  II.2 Bounds on the number of species

Before introducing axions in each sector, let us discuss the constraints on N.
These constraints arise from the bounds on the number of species $N_{s}$ due
to their interactions with gravity and via other forces. First, there exists a
gravitational bound on the number of particle species imposed by black hole
physics, which tells us that in a theory with $N_{s}$ particle species the
gravitational cutoff is not given by the Planck mass, but by the following
scale (usually called the species scale) Dvali:2007hz ; Dvali:2007wp ,

$M_{*}=\frac{M_{P}}{\sqrt{N_{s}}}\;.$ (2)

The scale $1/M_{*}$ gives the size of the smallest possible black hole.

Secondly, there exists another universal bound on species independent of
gravity. This bound says that the number of species $N_{s}$ inter-coupled by a
coupling of strength $\alpha(q)$ and the scale (momentum transfer) $q$ is
bounded by Dvali:2020wqi ; Dvali:2021jto

$N_{s}=\frac{1}{\alpha(q)}\;,$ (3)

where the coupling $\alpha$ must be evaluated at the scale $q$. As shown in
Dvali:2020wqi  this bound originates from species entropy since for a larger
number of species the unitarity would be violated due to multi-particle states
with exceedingly high entropy of species. In general, as shown there, no state
in the theory can have an entropy exceeding the value given by (3). The
coupling entering in the bound counts the subset of species that are
intercoupled. If the species are intercoupled selectively according to some
rules, the bound is modified accordingly. In each case, the guideline is the
entropy of species which leads to the conclusion that the coupling should not
exceed the bound given in (3). For explicit examples, see, Dvali:2020wqi . For
cosmological implications of the species entropy bound see Dvali:2023xfz .
Also, as explained there, from this perspective the black hole bound (2) on
species represents a particular case of species entropy bound (3), since the
quantity $(M_{*}/M_{P})^{2}$ represents a gravitational coupling evaluated at
the scale $p=M_{*}$.

In our analysis the minimal case would be a single extra axion ($N=2$ in
total) which corresponds to a single gauge group group in dark sector. In case
when the dark sector is represented by the exact SM copies, the minimal number
would be a single mirror copy. Possibility of dark matter in such a copy was
studied e.g., in Blinnikov:1982eh ; Kolb:1985bf . However, this value is
excluded if the dark copy is supposed to make up all the dark matter since in
this case, the dark matter halos would have collapsed into substructures. By
increasing the number of SM copies, the particle density in each sector is
reduced, leading to a suppressed dissipation efficiency and the lowest
phenomenologically viable value of $N$ being of $\mathcal{O}(10)$ Dvali:2007wp
; Dvali:2020wqi .

###  II.3 The Axion Mechanism

In order to eliminate the $\theta$-vacua let us now ad an axion WeinbergAxion
; Wilczek:1977pj  to each gauge sector. First, this can be done via a
traditional PQ formulation PQMechanism ; PQMechanism2 , which requires a
spontaneously broken global $U(1)_{PQ}$ symmetry. This means we introduce a
global U(1)PQ in each sector that is chiral, non-linearly realised and
anomalous with respect to the respective YM group. In the low-energy effective
theory the axion fields $a_{i}(x)$ are described as pseudo-scalars with the
anomalous couplings

$\mathcal{L}^{\rm
anomalous}_{i}\propto\left(\frac{a_{i}}{f^{a_{i}}}+\theta_{i}\right)G_{i}\tilde{G}_{i}\;,$
(4)

where the $G^{i}$ denote the YM field strength tensors, $\tilde{G}^{i}$ their
duals, and $f^{a_{i}}$ the axion decay constants. Below the confinement scale,
these couplings transform into potentials for the axions, which according to
the Vafa-Witten theorem have their vacua at $\langle
a_{i}\rangle=-f^{a_{i}}\theta_{i}$. Expanding around these vacua removes the
vacuum angles from the theory. Note that in more complex theories each YM
vacuum angle receives additional contributions, e.g. from the phases of the
quark matrices in the SM. We do not introduce an additional symbol for this
and keep the different contributions implicitly encoded in $\theta_{i}$.

The implementation of the U(1)PQ with the required properties into the SM is
usually achieved by the KSVZ and DFSZ invisible axion models DFSZ1 ; DFSZ2 ;
KSVZ1 ; KSVZ2 . The KSVZ model introduces an SM scalar singlet together with a
heavy quark, while the DFSZ model uses a singlet and an additional Higgs
doublet. In the former, the singlet $\Phi_{1}$ acquires a vacuum expectation
value $f^{a_{1}}$ via a standard Mexican hat potential,

$V_{1}(\Phi_{1})=\lambda_{1}\left(|\Phi_{1}|^{2}-(f^{a_{1}})^{2}\right)^{2}\;$
(5)

which generates the spontaneous breaking of the U(1)PQ. The axion is then
nothing else than the angular degree of freedom of the singlet, which we
canonically normalize as $\theta_{1}(x)\equiv a_{1}(x)/f^{a_{1}}$. In the SM
the introduction of the singlet is necessary to decouple the electroweak scale
from the PQ scale $f^{a_{1}}$, thus avoiding star cooling bounds and direct
detection until now. Here we focus on the minimal invisible axion models but
in principle, there exists a plethora of KSVZ- and DFSZ-type models, which
could be realized in our sector (see DiLuzio:2017pfr ; Plakkot:2021xyx  and
DFSZ  for catalogs of preferred models in each category).

The massless quark solution to the strong CP problem is in fact also a
realization of the PQ mechanism. In the presence of a massless quark, say the
up-quark in the SM, U(1)PQ symmetry coincides with the U(1)A symmetry of QCD
and the $\eta^{\prime}$ meson plays the role of the axion dvali2005threeform .
This implementation is not realized in the SM since all SM quarks are massive,
but it can be realized in a dark sector.

However, the axion mechanism can alternatively be realized via an intrinsic
gauge formulation, originally introduced in dvali2005threeform  and deeper
investigated indvali2022strong  (see also PhysRevD.105.085020  for a
consistency check). This formulation does not involve any global $U(1)_{PQ}$
symmetry. Instead, it is based on gauge redundancy of the respective QCD like
theory. Thus, in gauge formulation, the axion is introduced as a two-form
$B^{1}_{\mu\nu}$ that becomes the longitudinal mode of the QCD Chern-Simons
three-form $C^{1}_{\mu\nu\rho}$ (similar to the scalar that becomes the
longitudinal mode of the Maxwell-field in the Stückelberg formulation of the
Proca theory). This fixes the form of the theory, to the lowest order, as,

$\mathcal{L}_{\rm eff}=-\frac{1}{2\cdot
4!}(F^{1}_{\mu\nu\rho\sigma})^{2}+\frac{m^{2}}{2\cdot
3!}\left(C^{1}_{\mu\nu\rho}-\partial_{[\mu}B^{1}_{\nu\rho]}\right)^{2}\;,$ (6)

where $F^{1}_{\mu\nu\rho\sigma}\equiv\partial_{[\mu}C^{1}_{\nu\rho\sigma]}$
denotes the canonically normalized field strength tensor of the Chern-Simons
three-form. This puts the Chern-Simons theory into the Higgs phase,
effectively rendering the vacuum angle unphysical. The clear advantage of this
implementation is that it does not suffer from the quality problem, meaning
that $\theta_{1}=0$ can not be spoiled by any explicit breaking, since no
global symmetry is required. Correspondingly, the mechanism is stable under
any gauge invariant deformation of the theory. As discussed in dvali2022strong
, one of the consequences of this realization is that any measurement of a
neutron electric dipole moment would imply new physics, in contrast to the
ordinary PQ implementation where a non trivial $\theta_{1}$ could come from
any explicit breaking. Below $f^{a_{1}}$ the pseudo-scalar $a_{1}(x)$ and the
two-form axion $B^{1}_{\mu\nu}$ are dual to each other so that the physics is
vastly the same apart from the potential explicit breaking in ordinary PQ
realization. Above $f^{a_{1}}$, however, the formulations are very different
and the gauge formulation likely requires its UV-completion directly in
quantum gravity. This implies a PQ scale of $f^{a_{1}}\sim M_{\mathrm{P}}$ (or
$f^{a_{1}}\sim M_{*}$ in the presence of a large number of species).
Obviously, the gauge axion realization can be implemented for arbitrary hidden
sector gauge groups.

In the dark sectors, the concrete implementation of the U(1)PQ depends on the
model but in principle, there is more freedom due to a lack of bounds. For
exact SM copies, we could again use the KSVZ or DFSZ models or the two-form
implementation. For pure YM sectors, while we could again implement the axion
mechanism by using the KSVZ scenario or the gauge two-form implementation, the
massless quark solution is also viable.

In general early cosmologies of gauge and $U(1)_{PQ}$ axions can be very
different, in particular, due to fundamental differences in their UV
completions above the scale $f_{a}$. However, our initial conditions make
these potential differences unimportant. Hence, we shall remain independent of
any concrete realization of the axion mechanism and treat $f^{a_{i}}$ as free
parameters. However, if necessary, we will discuss consequences that are
specific to certain realizations.

###  II.4 Misalignment: Exact SM Copies

In the single axion case a very stringent bound on the axion scale results
from the misalignment mechanism being required not to produce more dark matter
than observed MisalignmentPRESKILL ; MisalignmentDINE ; MisalignmentAbbott .
An important input quantity in this bound is the initial misalignment angle
$\theta^{\rm ini}$, which can be defined as the value in the latest epoch when
the Hubble parameter crosses below the axion mass. Naturally, the most
stringent bound is obtained for $\theta^{\rm ini}\sim\mathcal{O}(1)$. However,
as shown in Dvali:1995ce  this value is highly sensitive to the details of the
previous evolution of the axion field (e.g., to the value of the QCD scale
during inflation) and can be made arbitrarily small via the early relaxation
mechanism. Due to this, in our analysis, we shall treat $\theta^{\rm ini}$ in
each gauge sector as independent free parameters. Under these circumstances,
we shall derive constraints on the parameter space for models with multiple
axions.

For a clearer presentation, let us begin with the case of $N$ exact SM copies,
where the equality between the sectors can be guaranteed by a discrete
symmetry Dvali:2007hz ; Dvali:2007wp . Various phenomenological aspects of
this framework were studied in FOOT199167 ; Dvali:2009ne . As already
explained Dvali:2007iv , despite the exact permutation symmetry of the
Lagrangian, the solution of the small $\theta$ puzzle requires the existence
of $N$ axions. The same is demanded by the S-matrix consistency of gravity
dvali2022strong . This necessitates the existence of $N$ axions. For
alternative compositions of the dark sectors we then only discuss the
differences in the calculation and give the final result.

In the standard axion picture, the U(1)PQ symmetry is spontaneously broken at
$T^{\rm PQ}\sim f_{a}$. In our framework, this takes place in every sector at
the temperatures $T^{\rm PQ}_{i}\sim f_{a_{i}}\equiv f_{a}$. The axions
$\theta_{i}(x)$ emerge as the corresponding Goldstone bosons and, as such,
have flat potentials. Using the standard notation we canonically normalize the
Goldstone-fields as $\theta_{i}(x)\equiv a_{i}(x)/f_{a}$.

The potential of each axion is flat only when the corresponding QCD is in the
Coulomb phase. In the confinement phase, each axion receives an effective
potential from non-perturbative effects of its corresponding QCD, which to
leading order in the semi-classical approximation can be accounted by
instantons. In a cosmological setting, these instantons are coupled to the
thermal bath, which causes the axion masses to depend on temperature. For
instance, in the dilute instanton gas approximation PhysRevD.17.2717 , the
potential takes the form QCDINstantonsFiniteTemp

$V_{i}(\theta_{i})=m_{a}^{2}(T_{i})f_{a}^{2}\big{(}1-\cos(\theta_{i})\big{)}\;,$
(7)

where the Temperature-dependent masses are given by

$m_{a}(T_{i})\equiv\frac{(\Lambda_{\rm
QCD}^{3}m_{u})^{\frac{1}{2}}}{f_{a}}\begin{cases}\beta\,\left(\frac{\Lambda_{\rm
QCD}}{T_{i}}\right)^{4}&:T_{i}>\Lambda_{\rm QCD}\;,\\\
1&:T_{i}\lesssim\Lambda_{\rm QCD}\;.\end{cases}$ (8)

Here, $m_{u}$ denotes the up-quark mass while $\beta$ encodes QCD and active
quark physics at the temperature $T_{i}$. For the SM, we roughly have
$\beta\sim 10^{-2}$ PhysRevD.17.2717 .

In an FLRW background with $R(t)$ denote the scale factor and $H(t)$ the
Hubble parameter, the equation of motion for each axion takes the form

$\ddot{\theta}_{i}+3H(t)\dot{\theta}_{i}-\frac{1}{R^{2}(t)}\Delta\theta_{i}+\frac{V_{i}^{\prime}(\theta_{i})}{f_{a}^{2}}=0\;.$
(9)

Let us use the potential (7) and make the following two simplifications.
First, only the leading order in $\theta_{i}$ is considered in the potential.
Secondly, only the zero mode of $\theta_{i}$ is taken into account. With these
simplifications, each equation of motion reduces to that of a damped harmonic
oscillator,

$\ddot{\theta}_{i}+3H(t)\dot{\theta}_{i}+m_{a}^{2}\big{(}T_{i}(t)\big{)}\theta_{i}=0\;.$
(10)

At $T_{i}>>\Lambda_{\rm QCD}$ the potential is flat so that
$\theta_{i}\approx\mathrm{const}$, i.e. the i-th axion is basically frozen
out. At $T_{i}\sim\Lambda_{\rm QCD}$ the corresponding axion potential gets
significant but it is not until the Hubble friction is overcome before the
axion starts performing coherent oscillations around the vacuum. We define
this critical time via

$m_{a}\big{(}T_{i}(t_{\rm osc})\big{)}=3H(t_{\rm osc})\;.$ (11)

From that moment on its equation of state no longer correspond to that of dark
energy but to that of non-relativistic matter, making the corresponding axion
contribute to the dark matter energy density.

For our model, this means that, due to the very low temperature in the dark
sectors, the axions there are essentially created with their zero-temperature
potential switched on. In contrast, our axion is created with a flat
potential. Therefore, since the zero-temperature mass is larger than the high-
temperature mass, the axion from our sector will start oscillating later than
the other axions. In order to quantify this, let us use our sector’s
temperature $T_{1}\equiv T$ as a clock instead of the cosmic time $t$. In
terms of $T$ the Hubble parameter during radiation domination is given by

$H(T)=\sqrt{\frac{\rho_{\rm
tot}}{3M_{\mathrm{P}}^{2}}}\sim\frac{T^{2}}{M_{\mathrm{P}}}\;.$ (12)

Here, $\rho_{\rm tot}$ denotes the total energy density, which is dominated by
our sector. Using (12) and the axion masses defined by (8), the condition (11)
results in the oscillations commencing when

$\displaystyle T_{{\rm osc},1}$ $\displaystyle\sim\left[\frac{\beta
M_{\mathrm{P}}\Lambda_{\rm
QCD}^{\frac{11}{2}}m_{u}^{\frac{1}{2}}}{3f_{a}}\right]^{\frac{1}{6}}$
$\displaystyle\sim 4\times
10^{-1}\,\mathrm{GeV}\left(\frac{10^{12}\,\mathrm{GeV}}{f_{a}}\right)^{\frac{1}{6}}\;,$
(13)

in our sector and

$\displaystyle T_{{\rm osc},i}$
$\displaystyle\sim\left[\frac{M_{\mathrm{P}}\Lambda_{\rm
QCD}^{\frac{3}{2}}m_{u}^{\frac{1}{2}}}{3f_{a}}\right]^{\frac{1}{2}}$
$\displaystyle\sim 2\times
10^{1}\,\mathrm{GeV}\left(\frac{10^{12}\,\mathrm{GeV}}{f_{a}}\right)^{\frac{1}{2}}\;,$
(14)

in the dark sectors with $i\neq 1$. Note that both moments in time are
expressed in terms of our sector’s temperature.

While we left $f_{a}$ as a free parameter there is a small caveat for
$f_{a}\gtrsim 6\times 10^{17}\,\mathrm{GeV}$. In this range our axion has also
reached its zero-temperature mass before overcoming the Hubble friction
Fox:2004kb , so that for these values of $f_{a}$ all axions start oscillating
at the temperature dictated by (14).

The initial energy density of each oscillation is

$\rho_{a_{i}}(T_{{\rm osc},i})=\frac{1}{2}f_{a}^{2}m_{a_{i}}^{2}(T_{{\rm
osc},i})\theta^{2}_{i}(T_{{\rm osc},i})\;,$ (15)

where $\theta_{i}(T_{{\rm osc},i})\equiv\theta_{i}^{\rm ini}$ is each sectors
initial misalignment angle. Essentially, we define the initial misalignment
angle $\theta_{i}^{\rm ini}$ as the value when the axion mass crosses below
the Hubble parameter in most recent history. However, as shown in Dvali:1995ce
this quantity is strongly sensitive to the pre-history and can be made
arbitrarily small due to early relaxation mechanisms. We thereby shall treat
$\theta_{i}^{\rm ini}$ as a free parameter.

Since Big Bang Nucleosynthesis took place during radiation domination, the
dark axions must not dominate the energy density of the universe at $T_{{\rm
osc},i}$. The energy density for relativistic degrees of freedom at that time
is given by $\rho_{\rm rad}(T_{{\rm osc},i})\sim g_{*}(T_{{\rm
osc},i})(T_{{\rm osc},i})^{4}$, where the effective number of relativistic
species is $g_{*}(T_{{\rm osc},i})\sim 10^{2}$. The requirement
$\sum_{i=2}^{N}\rho_{a_{i}}(T_{{\rm osc},i})\ll\rho_{\rm rad}(T_{{\rm
osc},i})$ results in the bound

$N\lesssim\frac{g_{*}(T_{{\rm osc},i})T_{{\rm osc},i}^{4}}{\Lambda_{\rm
QCD}^{3}m_{\mathrm{u}}(\theta_{i}^{\rm ini})^{2}}\sim
10^{12}\left(\frac{10^{12}\,\mathrm{GeV}}{f_{a}}\right)^{2}\left(\frac{1}{\theta_{i}^{\rm
ini}}\right)^{2}\;.$ (16)

As it turns out, a stronger bound on $N$ can be found at today’s Temperature.
Since the axions are decoupled, the number of zero modes per co-moving volume
is conserved as long as the changes in the mass are in the adiabatic regime.
Assuming this to be the case in every sector, today’s energy density per
sector is

$\rho_{a_{i}}(T_{\rm today})=\rho_{a_{i}}(T_{{\rm
osc},i})\frac{m_{a_{i}}(T_{\rm today})}{m_{a_{i}}(T_{{\rm
osc},i})}\left(\frac{T_{\rm today}}{T_{{\rm osc},i}}\right)^{3}\;$ (17)

We normalize with respect to the critical energy density $\rho_{\rm cr}\sim
M_{\mathrm{P}}^{2}H_{\rm today}^{2}$ to receive the corresponding axion
fraction $\Omega_{a_{i}}$ in today’s universe. By comparing $\Omega_{a_{i}}$
with the dark matter fraction $\Omega_{\rm DM}$ from the latest Planck mission
Aghanim:2018eyx , we find for our sector

$\frac{\Omega_{a_{1}}}{\Omega_{\rm DM}}\sim
0.54\left(\frac{\beta^{-\frac{1}{6}}m_{u}^{\frac{5}{12}}\Lambda_{\rm
QCD}^{\frac{7}{12}}}{10^{-2}\,\mathrm{GeV}}\right)\left(\frac{f_{a}}{10^{12}\,\mathrm{GeV}}\right)^{\frac{7}{6}}\left(\frac{\theta_{1}^{\rm
ini}}{1}\right)^{2}\;,$ (18)

while for each dark sector, i.e. $i\neq 1$, we get

$\frac{\Omega_{a_{i}}}{\Omega_{\rm DM}}\sim
0.01\left(\frac{m_{u}^{\frac{1}{4}}\Lambda_{\rm
QCD}^{\frac{3}{4}}}{10^{-2}\,\mathrm{GeV}}\right)\left(\frac{f_{a}}{10^{12}\,\mathrm{GeV}}\right)^{\frac{3}{2}}\left(\frac{\theta_{i}^{\rm
ini}}{1}\right)^{2}\;.$ (19)

In summary, the variations in reheating temperatures cause the axion
oscillations to commence at different times, resulting in varying axion
densities despite identical field theoretic parameters. While a single mirror
sector has an insignificant cosmological influence due to lower axion density
compared to our axion, larger numbers of mirror sectors result in the
accumulation of densities and their collective effect cannot be ignored. In
the following, we will examine the parameter space and explore the
consequences of this collective effect.

For our axion, a lower bound follows from astrophysics due to non-trivial
couplings to matter Raffelt:2006cw . These couplings lead to a contributions
in the cooling of stars that result in the bound $f_{a}\gtrsim
10^{9}\,\mathrm{GeV}$. Although other axions are not coupled to our sector
directly and are thus not contributing to the cooling in stars, the PQ scale
in all sectors must fulfill this bound due to the exact permutation symmetry
between the sectors.

Figure 1: The allowed parameter space with different initial misalignment
angles $\theta_{i}^{\rm ini}$ according to (21). The values on the thick lines
correspond to the case where the dark matter is entirely composed of axions.
The dashed line presents the species bound, meaning that along this line the
gravitational cutoff $M_{*}$ and $f_{a}$ coincide. Figure 2: The fraction of
our axion energy density with respect to the total axion energy density in the
case when the axions make up all the dark matter. For larger $f_{a}$ our axion
is dominating, while for lower $f_{a}$ the axions from dark sectors dominate.

While our axion density matter is not allowed to exceed the observed dark
matter density, leading to the known bound $f_{a}\lesssim
10^{12}\,\mathrm{GeV}$, the total axion density must not exceed it as well,
i.e.

$\sum_{i=1}^{N}\frac{\Omega_{a_{i}}}{\Omega_{\rm DM}}\lesssim 1\;.$ (20)

This results in the following inequality for the viable parameter space,

$\displaystyle N\lesssim$ $\displaystyle
10^{2}\left(\frac{10^{12}\,\mathrm{GeV}}{f_{a}}\right)^{3/2}\left(\frac{1}{\theta_{i}^{\rm
ini}}\right)^{2}$
$\displaystyle\times\left[1-0.54\left(\frac{f_{a}}{10^{12}\,\mathrm{GeV}}\right)^{\frac{7}{6}}\left(\frac{\theta_{1}^{\rm
ini}}{1}\right)^{2}\right]\;.$ (21)

We depict the viable regions in the $N$ \- $f_{a}$ \- plane for different
values of $\theta_{i}^{\rm ini}$ in Figure 1. The thick lines indicate the
combinations that result in the axions making up all the dark matter,
corresponding to the equality in (20). For example along the blue line in the
plot, the entire dark matter is composed of our axion when $f_{a}\sim
10^{12}\,\mathrm{GeV}$, whereas for $f_{a}\sim 10^{9}\,\mathrm{GeV}$
essentially the $N\sim 10^{6}$ axions from the dark sectors make up the dark
matter. A similar behaviour is observed for the other depicted values of
$\theta_{i}^{\rm ini}$.

We can quantify this by calculating the fraction of our axion along these
lines, which we depict in Figure 2. There we also see that for a fixed
$\theta_{i}^{\rm ini}$, lower $f_{a}$ values result in a dominant contribution
from dark sector axions, while larger $f_{a}$ values lead to a dominant
contribution from our axion.

Let us briefly discuss the range of the initial misalignment angles. The value
of $\theta_{i}^{\rm ini}$ depends on whether the PQ symmetry is broken during
or after inflation. When it is broken after inflation, for each sector there
appear all possible values of $\theta_{i}^{\rm ini}$ in today’s Hubble patch
so that one can average over the uniform distribution resulting in
$\theta_{i}^{\rm ini}=\pi/\sqrt{3}$ for all $i$. When the PQ symmetry is
broken during inflation, the initial misalignment angle is in principle an
arbitrary initial condition and $\theta_{i}^{\rm ini}\sim\mathcal{O}(1)$ seems
to be a reasonable choice in the absence of an explanation for special initial
conditions. However, such initial conditions appear naturally when QCD becomes
strong during inflation, yielding $\theta_{i}^{\rm ini}\ll 1$ Dvali:1995ce
(see also Takahashi:2018tdu ; Koutsangelas:2022lte  for concrete
implementations). If such a phase or a similar mechanism existed in each
sector, the parameter space would extend to $f_{a}\gg 10^{12}\,\mathrm{GeV}$
and $N\gg 10^{6}$. Note that for a given $N$, the PQ scale is bounded from
above by (LABEL:Eq:SpeciesBound), which we depict via a red dashed line in
Figure 1.

The existence of such a phase in every sector would allow for $N>10^{6}$ but
this poses the question if this is enough to solve the Hierarchy problem. As
mentioned in the previous section, the existence of $N\sim 10^{32}$ copies of
the SM would provide a solution to the hierarchy problem. As can be read of
from Figure 1 or equivalently from (21), this large number requires $f_{a}\ll
10^{9}\,\mathrm{GeV}$ together with an early phase of strong QCD. This is in
conflict with the lower bound on the PQ scale from astrophysical
considerations in our sector, i.e. $f_{a}\gtrsim 10^{9}\,\mathrm{GeV}$.
Seemingly this excludes a large number of exact SM copies as a solution for
the hierarchy problem. However, there are proposals such as the clockwork
mechanism Kaplan:2015fuy  that allow the separation of the PQ scale from the
suppression factor of the couplings. This would render the part of the
parameter space that solves the hierarchy problem viable.

###  II.5 Misalignment: Pure YM Dark Sector

Finally, let us discuss the case where each dark sector is based on a pure YM
group SU${}_{i}(N_{C})$. The crucial difference compared to the case with
exact SM copies would be the absence of light quarks. Without light quarks,
chiral perturbation techniques can no longer be used to calculate the zero-
temperature mass of the axion. Instead, the mass would have to be calculated
by using alternatives such as large $N_{c}$ methods. We use as an
approximation for the axion mass the extrapolated result of the dilute
instanton gas at finite temperatures QCDINstantonsFiniteTemp ,

$m_{a_{i}}(T^{i})\equiv\frac{(\Lambda^{i}_{\rm
conf})^{2}}{f_{a_{i}}}\begin{cases}\left(\frac{\Lambda^{i}_{\rm
conf}}{T^{i}}\right)^{4}&:T^{i}>\Lambda^{i}_{\rm conf}\;,\\\
1&:T^{i}\lesssim\Lambda^{i}_{\rm conf}\;,\end{cases}$ (22)

where $\Lambda_{{\rm conf},i}$ denotes the different confinement scales of the
dark sectors and the factor of $\beta$ was dropped due to the absence of light
quarks. In contrast to the case of exact SM copies, the $\Lambda_{{\rm
conf},i}$ are free parameters that in principle can be even smaller than the
dark sector temperatures. Hence, the dark sector axion can in principle be
created with an essentially flat potential, although this requires quite low
confining temperatures. In the following, we will only consider the case where
the dark axions are created with their zero-temperature mass.

Figure 3: The allowed region for a single dark, pure YM sector with
confinement scale $\Lambda_{\rm conf}$, axion scale $f_{a_{2}}$, and no
intersector interactions. Minimality favors $\theta^{\rm ini}_{2}\sim 1$,
since in this scenario there is no need for additional physics that results in
a small misalignment angle. The red line represents the perturbative unitarity
bound $\Lambda_{\rm conf}\lesssim f_{a_{2}}$.

From there the calculation is vastly the same as in the exact SM case. The
oscillations in the dark sectors commence at

$\displaystyle T_{{\rm osc},i}\sim\Lambda^{i}_{\rm
conf}\left(\frac{M_{\mathrm{P}}}{3f_{a_{i}}}\right)^{\frac{1}{2}}\;,$

which results in the final dark sector density of

$\frac{\Omega_{a_{i}}}{\Omega_{\rm DM}}\sim
10^{-18}\left(\frac{\Lambda^{i}_{\rm
conf}}{\,\mathrm{GeV}}\right)\left(\frac{\theta_{i}^{\rm
ini}}{1}\right)^{2}\left(\frac{f_{a_{i}}}{\,\mathrm{GeV}}\right)^{\frac{3}{2}}\;.$
(23)

Let us discuss the parameter space for the particularly interesting case of a
single, pure YM sector, i.e. $N=2$ with $i=2$ being the dark sector. To
simplify the notation, we define $\Lambda^{i=2}_{\rm conf}\equiv\Lambda_{\rm
conf}$. For the dark sector axion to not result in an abundance of dark
matter, we find

$\Lambda_{\rm conf}\lesssim
10^{18}\,\mathrm{GeV}\left(\frac{1}{\theta_{2}^{\rm
ini}}\right)^{2}\left(\frac{\,\mathrm{GeV}}{f_{a_{2}}}\right)^{\frac{3}{2}}\;,$
(24)

where the equality is valid when the dark axion makes up all the dark matter.
Furthermore, to ensure perturbative unitarity in the dark sector, we must have

$\Lambda_{\rm conf}\lesssim f_{a_{2}}\;.$ (25)

This follows from the requirement of the axion’s quartic self-interaction to
be in the weak coupling regime. We depict the viable region in Figure 3. If
the dark sector axion is supposed to make up the dark matter, we see that for
$\Lambda_{\rm conf}$ to lie below the Planck scale, the dark sector PQ scale
must fulfill $f_{a_{2}}\gtrsim 1\,\mathrm{GeV}$, which is pushed higher for
smaller values of $\theta_{2}^{\rm ini}$. Equivalently, for $f_{a_{2}}$ to lie
below the Planck scale, the dark confinement scale must fulfill $\Lambda_{\rm
conf}\gtrsim 1\,\mathrm{eV}$. This minimum value is also pushed higher for
smaller values of $\theta_{2}^{\rm ini}$.

It should be mentioned that the region with $\theta_{2}^{\rm ini}\sim 1$ is
favored by minimality in the dark sector, as tiny misalignment angels require
new physics Dvali:1995ce ; Koutsangelas:2022lte . Without a discrete symmetry
that connects $\Lambda_{\rm conf}$ and $f_{a_{2}}$ to the limited analogs in
our sector, there is simply no need to introduce new physics to achieve tiny
misalignment angles. Intriguingly, in the two-form realization, for which the
value of $f_{a_{2}}\sim M_{\mathrm{P}}$ is strongly favored, this minimality-
argument results in the dark confinement scale being $\Lambda_{\rm conf}\sim
1\,\mathrm{eV}$.

###  II.6 Isocurvature Perturbations

When one of the PQ symmetries is broken during (or before) inflation and never
restored afterwards, the corresponding axion field is subject to quantum
perturbations. Let us write the axion fields as
$\theta_{i}=\langle\theta_{i}\rangle+\delta\theta_{i}$, where
$\langle\theta_{i}\rangle=\theta_{i}^{\rm ini}$ and $\delta\theta_{i}$ denotes
the quantum perturbations. These fulfill $\langle\delta\theta_{i}\rangle=0$
and have a standard deviation of

$\sigma_{\theta_{i}}\sim\sqrt{\langle\delta\theta_{i}^{2}\rangle}\sim\frac{H_{\mathrm{I}}}{2\pi
f_{a_{i}}}\;.$ (26)

Furthermore, if none of the axions is responsible for inflation, their
fluctuations will not be of adiabatic but of isocurvature-type. Since these
lead to a unique imprint in the temperature and polarisation fluctuations of
the CMB, they give rise to a constraint on the axion’s parameter space.

One could worry that in the presence of several axions, the amount of
isocurvature perturbations is significantly enhanced, thus making the existing
constraint much more severe. Let us show that this is not necessarily the
case, even though the contributions are additive.

Assuming the perturbations to be normal distributed in the regime of small
$\theta_{i}$, where anharmonic corrections of the potentials can be ignored,
the collective amplitude of the axions isocurvature fluctuations is given by
Kobayashi:2013nva

$\displaystyle\Delta_{a}(k_{0})$ $\displaystyle=\frac{\delta\Omega_{\rm
DM}}{\Omega_{\rm DM}}=\frac{\sum_{i}\Omega_{a_{i}}}{\Omega_{\rm
DM}}\frac{\delta\ln\Omega_{a_{i}}}{\delta\theta_{i}^{\rm
ini}}\sigma_{\theta_{i}}$
$\displaystyle=\sum_{i}\frac{\Omega_{a_{i}}}{\Omega_{\rm
DM}}\frac{H_{\mathrm{I}}}{\pi\theta_{i}^{\rm ini}f_{a_{i}}}\;,$ (27)

where $\theta_{i}^{\rm ini}\ll\sigma_{\theta_{i}}$ was used. The latest
experimental bound on uncorrelated isocurvature perturbations by Planck is
Aghanim:2018eyx

$\beta(k_{0})\equiv\frac{\Delta^{2}_{a}(k_{0})}{\Delta^{2}_{\mathcal{R}}(k_{0})+\Delta^{2}_{a}(k_{0})}<0.038\quad{\rm
at\ 95\%\ CL}\;,$ (28)

where $k_{0}=0.050{\rm Mpc}^{-1}$. This can be translated to a constrain on
$H_{\mathrm{I}}$,

$H_{\mathrm{I}}\lesssim 10^{7}\,\mathrm{GeV}\sum_{i}\frac{\Omega_{\rm
DM}}{\Omega_{a_{i}}}\left(\frac{f_{a_{i}}}{10^{12}\,\mathrm{GeV}}\right)\left(\frac{\theta_{i}^{\rm
ini}}{1}\right)\;.$ (29)

For both models at hand, namely a single dark YM sector and $N$ exact SM
copies, the generated isocurvature perturbations are either dominated by our
sector or the dark sector(s).

For the former, if the dark matter is dominated by the axion from the dark
sector, we can express the bound on $H_{\mathrm{I}}$ in terms of the dark
confinement scale by using the equality of (24) to eliminate $\theta_{2}^{\rm
ini}$. This results in

$H_{\mathrm{I}}\lesssim
10^{7}\,\mathrm{GeV}\left(\frac{f_{a_{i}}}{10^{12}\,\mathrm{GeV}}\right)^{\frac{1}{4}}\left(\frac{1\,\mathrm{GeV}}{\Lambda_{\rm
conf}}\right)^{\frac{1}{2}}\;.$ (30)

In the case of a two-form axion from such a dark sector, i.e. with
$f_{a_{2}}\sim M_{\mathrm{P}}$ and the minimal possible dark confinement scale
$\Lambda_{\rm conf}\sim\,\mathrm{eV}$, the bound reduces to
$H_{\mathrm{I}}\lesssim 10^{13}\,\mathrm{GeV}$, meaning that such an axion is
essentially not plagued by an isocurvature constraint.

For $N$ exact SM copies, when the dark matter is collectively composed by the
axions from the copies, i.e. $\sum_{i}\Omega_{a_{i}}\sim\Omega_{\rm DM}$,
using (21) to eliminate the initial misalignment angles in (29) yields

$H_{\mathrm{I}}\lesssim
10^{7}\,\mathrm{GeV}\left(\frac{f_{a_{i}}}{10^{12}\,\mathrm{GeV}}\right)^{\frac{1}{4}}\left(\frac{10}{\sqrt{N}}\right)\;.$
(31)

We observe that for moderate values of $N\sim 10^{2}$ the bound remains
relatively unaffected. On the other hand, for values required to solve the
Hierarchy problem, i.e. $N\sim 10^{32}$ and $f_{a}\sim\,\mathrm{TeV}$, the
bound naively tightens to $H_{\mathrm{I}}\lesssim 10^{-10}\,\mathrm{GeV}$.
However, in this case, the dark confinement scale and the axion mass exceed
the inflationary Gibbons-Hawking temperature $T_{\mathrm{I}}\sim
H_{\mathrm{I}}$, meaning that the dark YM sector becomes strongly coupled
during inflation. Consequently, the axion develops a substantial mass so that
the requirements for the development of isocurvature perturbations are no
longer given. Of course, this requires avoiding the astrophysical bounds on
$f_{a}$ by a mechanism such as clockworking (see discussion at the end of
Section II.4).

###  II.7 $N$-Boson Stars

Bosonic dark matter is known to be able to form dense clumps, which in the
case of the axion are called axion stars. The expression “star” is used to
denote an object sustained by hydrostatic equilibrium, regardless of its
emission of light. Such objects as well as other variations of compact dark
matter structures have been long studied and provide an interesting
perspective on the nature of dark matter TKACHEV1991289  (see Zhang:2018slz
for a review).

For an axion arising from a dark YM sector, the physics of axion stars are
essentially equal to those of our sector, merely the values of the axion mass
and decay constant change. However, the existence of multiple dark sectors may
lead to structures that are collectively formed from the axions of the
different sectors (see CompactDarkMatterObjectsViaN  for a realization of this
idea with fermions called $N$-MACHOS). These structures would be much lighter
and smaller compared to the single-sector case, which can be intuitively
understood as follows. When particles from $N$ sectors are placed into a
volume, the system is subject to gravity and pressure. Gravity attempts to
collapse the system, while pressure from each sector resists the collapse.
Eventually, the combined pressure can counteract gravity, resulting in an
equilibrium configuration. However, because each particle can only interact
with a small portion of the total particles, the overall pressure is much
lower than in a single sector case. At the same time, gravity affects all
particles equally. As a result, equilibrium can only be achieved with smaller
masses and radii than in typical single-sector models. In the following, we
argue that this phenomenon can not take place for axions.

The analysis performed in CompactDarkMatterObjectsViaN  was based on the
Tolman-Oppenheimer-Volkoff (TOV) equation, which is a relativistic equation
that characterizes the hydrostatic equilibrium structure of a self-
gravitating, spherically symmetric object made of a perfect fluid. Unlike
perfect fluids, the scalar field of axion stars has a non-trivial self-
interaction potential that cannot be described by a simple equation of state.
As a result, the equilibrium structure of axion stars cannot be determined
solely by the TOV equation. Instead, it requires the solution of the Klein-
Gordon equation that describes the scalar field, which must be solved
simultaneously with the Einstein field equations.

Although a comprehensive analysis would certainly be of interest, here, we
limit ourselves to so called dilute axion stars. To identify such
configurations, we follow the standard approach (see for instance
Zhang:2018slz ) and employ a set of approximations to simplify the relevant
equations while still maintaining a high level of accuracy. We aim to maintain
a maximal level of generality in our calculations to allow for the application
to scalars with various potentials.

To begin, we can utilize non-relativistic effective field theory since the
axions produced through misalignment are non-relativistic. In this framework,
it is convenient to introduce the complex field $\psi_{i}(x)$ for each axion
through

$a_{i}(x)=\frac{1}{\sqrt{2m_{a_{i}}}}\left[\psi_{i}(x)\mathrm{e}^{im_{a_{i}}t}+\psi_{i}^{*}(x)\mathrm{e}^{-im_{a_{i}}t}\right]\;.$
(32)

By taking the non-relativistic limit of each Klein-Gordon equation and
expressing the axion field in terms of $\psi_{i}$, we obtain the time-
dependent Gross-Pitaevskii equation. To furthermore separate the time
dependence, we make the standard separation ansatz
$\psi_{i}(\vec{x},t)=A_{i}(\vec{x})\mathrm{e}^{-i\omega_{i}t}$. This results
in the time-independent Groß-Pitaevskii equation,

$m_{a_{i}}\left(\Phi+V_{i}^{\prime}(|\psi_{i}|^{2})\right)-\frac{1}{2m_{a_{i}}}\frac{\Delta\psi_{i}}{\psi_{i}}=E_{i}\;.$
(33)

Next, the cosmic axion condensate has a relatively low mean mass density, so
we can use Newtonian gravity described by the Poisson equation,

$\Delta\Phi=\frac{\sum_{i}m_{a_{i}}|\psi_{i}|^{2}}{M_{\mathrm{P}}^{2}}\;.$
(34)

Crucially, the sum appears because gravity couples to all the axions. Within
these approximations, the Groß-Pitaevskii-Poisson equations provide a
simplified set of equations that can accurately describe dilute axion stars.
We prefer to rewrite this system in the form of hydrodynamic equations. To do
so, we first take the Laplacian of the Gross-Pitaevskii equations and plug in
the Poisson equation. Secondly, we perform the Madelung transformation to
express everything in terms of the (pseudo) density
$\rho_{i}=m_{a_{i}}|\psi_{i}|^{2}$ and the pressure, which is related to the
potential via $\nabla V_{i}^{\prime}(\rho_{i})=\nabla
p_{i}(\rho_{i})/\rho_{i}$. This results in the fundamental equation of
hydrostatic equilibrium with quantum effects,

$-\nabla\cdot\left(\frac{\nabla
p_{i}}{\rho_{i}}\right)+\frac{1}{2m_{a_{i}}^{2}}\Delta\left(\frac{\Delta\sqrt{\rho_{i}}}{\sqrt{\rho_{i}}}\right)=\frac{\rho_{\rm
tot}}{M_{\mathrm{P}}^{2}}\;.$ (35)

The first term accounts for the hydrodynamic pressure, which can either be
attractive or repulsive, while the second term characterizes the repulsion
from quantum pressure due to the Heisenberg uncertainty principle. On the
right hand side the total density, $\rho_{\rm tot}=\sum_{i}\rho_{i}$, appears
because gravity couples to all axions.

The last approximation regards the (pseudo) density $\rho$, which we take to
be small compared to the cosmic condensate density
$m_{a_{i}}^{2}f_{a_{i}}^{2}$. In this regime the effective potential
$V_{i}(\rho_{i})$ is dominated by the leading term in its power series. For
the instantonic potential from (7), the leading order interaction is the
quartic term of $a_{i}$, so in terms of $\rho_{i}$ the potential can be
approximated as

$V_{i}(\rho_{i})\sim-\frac{1}{m_{a_{i}}^{2}f_{a_{i}}^{2}}\rho_{i}^{2}\;.$ (36)

This is equivalent to the negative polytropic pressure

$p_{i}\sim-\frac{1}{m_{a_{i}}^{2}f_{a_{i}}^{2}}\rho_{i}^{2}\equiv
K_{i}\,\rho_{i}^{\gamma_{i}}\;,$ (37)

with polytropic constant $K_{i}\sim-1/(m_{a_{i}}^{2}f_{a_{i}}^{2})$ and
polytropic index $\gamma_{i}=2$. With such an equation of state, (35) becomes

$-\nabla\left(\gamma_{i}K_{i}\rho_{i}^{\gamma_{i}-2}\nabla\rho_{i}\right)+\frac{1}{2m_{a_{i}}^{2}}\Delta\left(\frac{\Delta\sqrt{\rho_{i}}}{\sqrt{\rho_{i}}}\right)=\frac{\rho_{\rm
tot}}{M_{\mathrm{P}}^{2}}\;.$ (38)

For the sake of illustration, let us assume that all densities are equal,

$\rho_{i}=\rho_{\rm tot}/N\;,$ (39)

as collective effects are most pronounced when all sectors have equal
densities CompactDarkMatterObjectsViaN . In this way, we can express (38) in
terms of the total density,

$-\nabla\left(\frac{\gamma_{i}K_{i}\rho_{\rm
tot}^{\gamma_{i}-2}}{N^{\gamma_{i}-1}}\nabla\rho_{\rm
tot}\right)+\frac{1}{2m_{a_{i}}^{2}}\Delta\left(\frac{\Delta\sqrt{\rho_{\rm
tot}}}{\sqrt{\rho_{\rm tot}}}\right)=\frac{\rho_{\rm
tot}}{M_{\mathrm{P}}^{2}}\;.$ (40)

We observe that in the presence of $N$ axions the pressure from short-range
interactions is suppressed by $N^{\gamma_{i}-1}$ while the quantum pressure is
unaffected.

The influence of $N$ on the solution can be understood by considering the
following limiting cases.

  * •

Neglecting the hydrostatic pressure, the equilibrium is between the
gravitational attraction and the repulsion by the quantum pressure, which are
both independent of $N$. Thus, the final configuration is unchanged from the
single axion case.

  * •

Neglecting the quantum pressure, the equilibrium is between the gravitational
attraction and the hydrostatic pressure (which must be repulsive), which is
now strongly suppressed by $N^{\gamma_{i}-1}$. From here on the analysis is
exactly the same as in CompactDarkMatterObjectsViaN , so that the radius and
mass of the final configuration would be suppressed by $1/\sqrt{N}$.

  * •

Neglecting the gravitational attraction, the equilibrium is between the
hydrostatic pressure (which must be attractive) and the repulsion by quantum
pressure. This would again result in an altered mass spectrum, however, such
equilibria are always unstable Chavanis_2011 .

We can see that stable collective configurations with a suppressed mass
spectrum can only exist for repulsive interactions, which is not the case for
the axion. Therefore, dilute axion stars would consist of several axions but
would not differ in mass or radius compared to the single sector case.
However, we want to stress that they would exist for scalars with a repulsive
interaction (see for instance Berezhiani:2021rjs ; Berezhiani:2022buv )

##  III Kinetic Mixing between Axions

###  III.1 Kinetic Mixing as Intersector Interaction

In this section, we discuss non-gravitational intersector interactions by
restoring $\mathcal{L}_{\rm mix}$ in (1). For now, we again put the focus on
$N$ exact SM copies. Without the additional axions, the possible
renormalizable interactions of this model, which are compatible with gauge-,
Lorentz-, and the underlying discrete symmetry, are photon kinetic mixing, a
Higgs portal coupling, and neutrino mass mixing. At the non-renormalizable
level, additional interactions such as neutron oscillations become possible as
well. The effects of these interactions have been vastly discussed in the
literature (see Foot:2014mia  for a review). As our focus is not on these
interactions, we will not elaborate further on their influence. Instead, we
focus on axion kinetic mixing, which in the model under consideration is
described by

$\mathcal{L}_{\rm mix}=\epsilon\sum_{i\neq
j}\partial_{\mu}a^{i}\partial^{\mu}a^{j}\;,$ (41)

where $\epsilon$ parametrizes the kinetic mixing strength. An important
implication of kinetic mixing between our visible and hidden sectors is that
it can result in thermal equilibrium between the sectors. In order to prevent
inconsistencies with nucleosynthesis or an excessive amount of dark matter
particles, the thermal equilibrium must not occur prior to BBN, usually
leading to a constraint on the kinetic mixing strength. The axions that are
produced by misalignment are non-thermal. Thus, they can not transfer heat
between the sectors via axion-axion interactions to thermalize any sectors.
However, this is not true for other types of interactions, potentially leading
to a constraint on $\epsilon$ in order to avoid reheating of the dark sectors
prior to nucleosynthesis. Such an analysis, which requires further
assumptions, goes beyond the scope of this paper and we simply regard
$\epsilon$ as a (phenomenologically) free parameter in the following.

The way in which axion kinetic mixing can emerge in the low energy effective
theory depends on the specific implementation of the axion. In the standard
implementation of the axion via the KSVZ model a singlet Higgs is added to
make the axion invisible. Denoting this singlet in each sector as $\Phi_{i}$
allows for the following dimension six Higgs portal couplings between the
sectors

$\mathcal{L}_{\rm mix}^{\rm UV}=\frac{1}{M^{2}}\sum_{i\neq
j}(\Phi_{i}^{\dagger}\partial_{\mu}\Phi_{i})(\Phi_{j}\partial^{\mu}\Phi_{j}^{\dagger})+{\rm
h.c.}\;,$ (42)

where $M$ is some cut-off scale. The singlets acquire the VEVs
$f_{a_{i}}\equiv f_{a}$ via a proper scalar potential that spontaneously
breaks the U(1)PQ symmetries. This results in the kinetic mixing described in
(41) with

$\epsilon\equiv\frac{f_{a}^{2}}{2M^{2}}\;.$ (43)

A diagram coming from such a operator is shown in Figure 4. Similar operators
that result in axion kinetic mixing are also possible in DFSZ-type models,
where extra Higgs doublets are required in addition to the singlet.

While we neglect the phenomenological constraint on $\epsilon$ arising from
the necessary avoidance of the dark sector thermalization before BBN, it is
still subject to consistency limitations due to unitarity. The operators in
(42) already at tree-level put a unitarity constraint:

$\epsilon\lesssim\frac{1}{\sqrt{N}}\;.$ (44)

However, the bound on $\epsilon$ is much more severe. Since it describes an
intersector interaction coupling at momentum transfer $q\sim f_{a}$ among all
species, it has to fulfill the bound (3) of such a coupling imposed by the
entropy of species Dvali:2020wqi . The bound (3) gives:

$\epsilon\lesssim\frac{1}{N}\;.$ (45)

In the alternative two-form implementation of the axion, axion kinetic mixing
appears via non-diagonal mass terms of three-form action,

$\mathcal{L}_{\rm mix}\sim\sum_{i\neq
j}\left(C^{i}_{\mu\nu\rho}-\partial_{[\mu}B^{i}_{\nu\rho]}\right)\left(C^{j}_{\mu\nu\rho}-\partial_{[\mu}B^{j}_{\nu\rho]})\right)\;.$
(46)

Mixing of this kind is not forbidden by gauge symmetries and can originate
from virtual black hole exchange Dvali:2007iv , but there may be other
origins. The emergence of such operators is for $i\neq j$ only possible if
they are suppressed by powers of $M_{\mathrm{P}}^{-1}$. The reason for this
suppression is that micro black holes cannot be universally coupled, which
necessitates gravitational suppression of inter-sector transitions at the
fundamental level.

Figure 4: Loop induced by the operators in (42). Each vertex contributes the
effective coupling $\Delta p^{2}/M^{2}$, but for maximum momentum transfers of
order $f_{a}$ the effective coupling becomes $\epsilon$.

In the remainder of this section, we will study the consequences of the axion
kinetic mixings in the pseudo-scalar formulation. We will stay agnostic of
their origin, but use the unitarity bound described by (45).

###  III.2 Modification of Axion Physics for SM Copies

The presence of the extra term in the Lagrangian necessitates a change of
basis to express it in a canonical form. As a result, there is a mismatch
between the “sector basis”, which corresponds to the labels of the species,
and the “canonical kinetic basis”, in which the propagator is in canonical
form. So the task is to find the relation between these two bases.

Before discussing this change of basis, it is worth noting that a similar
situation arises in the SM’s neutrino sector. There, the interaction terms are
diagonal in the flavor basis. But in this basis the mass matrix of neutrinos
is off-diagonal. Therefore, there is a mismatch between the flavor and the
mass basis, which leads to phenomenons like neutrino oscillations. In our
case, the mismatch appears between the canonical kinetic basis and the sector
basis.

The first step is to express the kinetic part in (41) as

$\mathcal{L}\supset\begin{pmatrix}\partial_{\mu}a_{1}\\\ \vdots\\\ \vdots\\\
\partial_{\mu}a_{N}\end{pmatrix}^{T}\begin{pmatrix}1&&\epsilon&&\dots&&\epsilon\\\
\epsilon&&1&&\ddots&&\vdots\\\ \vdots&&\ddots&&\ddots&&\epsilon\\\
\epsilon&&\dots&&\epsilon&&1\end{pmatrix}\begin{pmatrix}\partial_{\mu}a_{1}\\\
\vdots\\\ \vdots\\\ \partial_{\mu}a_{N}\end{pmatrix}\;.$ (47)

This matrix appears due to the permutation symmetry among the copies also in
other cases, which have been studied in Dvali:2009ne . We can rewrite the
matrix, which we shall call $K$ from now on, in the following way,

$K=\begin{pmatrix}1-\epsilon&0&\dots&0\\\ 0&1-\epsilon&\ddots&\vdots\\\
\vdots&\ddots&\ddots&0\\\
0&\dots&0&1-\epsilon\end{pmatrix}+\epsilon\begin{pmatrix}1&\dots&1\\\
\vdots&\ddots&\vdots\\\ 1&\dots&1\end{pmatrix}\;.$ (48)

The problem reduces then to the diagonalization of a matrix of just ones. The
transformation matrix $S$ that diagonalizes $K$ is

$S=\begin{pmatrix}1&1&1&\dots&1\\\ 1&-1&0&\dots&0\\\ 1&0&-1&\ddots&\vdots\\\
\vdots&\vdots&\ddots&\ddots&\vdots\\\ 1&0&\cdots&0&-1\end{pmatrix}\;.$ (49)

The first row gives rise to the eigenstate of the form,

$\tilde{a}_{L}=\frac{1}{\sqrt{N}}a_{1}+\sqrt{\frac{N-1}{N}}a_{h}\;,$ (50)

where we have introduced the notation

$a_{h}=\frac{1}{\sqrt{N-1}}\sum_{i=2}a_{i}\;,$ (51)

due to later convenience. The eigenstate $\tilde{a}_{H}$ corresponds to the
eigenvalue of $1+(N-1)\epsilon$. Because the matrix $S$ is not a unitary
matrix we still have to find a convenient basis.

From all other rows, we see that we have $N-1$ degenerate eigenstates $v_{i}$
of the eigenvalue $1-\epsilon$. These $v_{i}$´s are the columns of $S$. Due to
this degeneracy, we can reduce the $N\times N$ problem to a $2\times 2$
problem by defining a superposition of these degenerate eigenstates. This
looks like

$\tilde{a}_{H}=\frac{1}{N-1}\sum_{i}v_{i}=a_{1}-\frac{1}{\sqrt{N-1}}a_{h}\;,$
(52)

and again has the eigenvalue $1-\epsilon$. After normalization this becomes

$\tilde{a}_{H}=\sqrt{\frac{N-1}{N}}a_{1}-\frac{1}{\sqrt{N}}a_{h}\;.$ (53)

We see that $\tilde{a}_{H}$ is a collective expression made out of all former
$a_{i}$´s. Using expressions (53) and (50), we find $a_{1}$ expressed in the
canonical kinetic basis

$a_{1}=\sqrt{\frac{N-1}{N}}\tilde{a}_{H}+\frac{1}{\sqrt{N}}\tilde{a}_{L}\;.$
(54)

Note that the choice of what we label as species one is totally arbitrary and
therefore holds for every $i=1,...,N$. This results from having the same
parameter $\epsilon$ among all sectors. The $2\times 2$ matrix connecting
$a_{1}$ and $a_{h}$ with $\tilde{a}_{H}$ and $\tilde{a}_{L}$ is now a unitary
matrix.

Expressing (47) in terms of $\tilde{a}_{L}$ and $\tilde{a}_{H}$, there is no
more mixing but each term is multiplied by the corresponding eigenvalue. In
order to have canonical kinetic terms, the states $\tilde{a}_{L}$ and
$\tilde{a}_{H}$ need to be redefined by

$\tilde{a}_{H}\rightarrow\frac{1}{\sqrt{1-\epsilon}}\tilde{a}_{H}\;,\quad\tilde{a}_{L}\rightarrow\frac{1}{\sqrt{1+(N-1)\epsilon}}\tilde{a}_{L}\;.$
(55)

One effect of this redefinition is that the mass part in the Lagrangian is
changed to

$\mathcal{L}_{\rm mass}\sim\begin{pmatrix}\tilde{a}_{H}\\\
\tilde{a}_{L}\end{pmatrix}^{T}\begin{pmatrix}\frac{m_{a}^{2}}{1-\epsilon}&&0\\\
0&&\frac{m_{a}^{2}}{1+(N-1)\epsilon}\end{pmatrix}\begin{pmatrix}\tilde{a}_{H}\\\
\tilde{a}_{L}\end{pmatrix}\;.$ (56)

with $m_{a}$ being the mass induced by the PQ-mechanism. One sees that kinetic
mixing leads to a splitting of the masses of the axions. $N-1$ axions are
degenerated and one light axion whose mass is suppressed by the number of
copies. In other words, kinetic mixing of many copies of the axion leads to
two different detectable axion states which have different masses. The
relation between these masses is

$\frac{m_{L}}{m_{H}}=\sqrt{\frac{1-\epsilon}{1+(N-1)\epsilon}}\sim\frac{1}{\sqrt{2}}\;,$
(57)

where for the approximation we assumed natural values of $\epsilon$ with
respect to the unitarity bound, i.e. $\epsilon\sim 1/N$, and $N\gg 1$.

Figure 5: Viable axion-photon couplings for the light axion state (yellow
band), which arises from equal kinetic mixing of axions from $N$ sectors.
Current bounds and predicted sensitivities of future experiments are depicted
in different colors. The parameter space for $N\lesssim 10$ will be probed by
future experiments. The yellow band is not to be confused with the bands from
DFSZ , which include the set of DFSZ-type or KSVZ-type models. The two stars
represent a situation when future experiments measure one axion (the upper
star) and where we expect the second axion (the lower star). This plot was
created with the help of the software AxionLimits .

With the relation between the species basis and the canonical kinetic basis,
we can turn to the phenomenological implications of this model. Due to the
Goldstone nature of the axion, its lowest order couplings to fermions and
gauge bosons have the generic form

$g_{ao}a_{1}\mathcal{O}\;,$ (58)

where the whole UV dependency is encoded in $g_{ao}$. After the
diagonalization and field redefinition of (55), the coupling becomes,

$g_{ao}\left(\sqrt{\frac{N-1}{N}}\frac{1}{\sqrt{1-\epsilon}}\tilde{a}_{H}+\frac{1}{\sqrt{N}}\frac{1}{\sqrt{1+(N-1)\epsilon}}\tilde{a}_{L}\right)\mathcal{O}\;.$
(59)

We observe that by having $N$ kinetically mixed axions, our sector couples to
$N$ axions instead of one. These $N$ axions come in two categories: $N-1$
axions ($\tilde{a}_{H}$) that behave exactly the same, and one special axion
($\tilde{a}_{L}$). The single axion coupling $g_{ao}$ gets modified by one of
the following factors,

$f_{H}(N,\epsilon)\equiv\sqrt{\frac{N-1}{N}}\frac{1}{\sqrt{1-\epsilon}}\sim
1\;,\\\ $ (60)
$f_{L}(N,\epsilon)\equiv\frac{1}{\sqrt{N}}\frac{1}{\sqrt{1+(N-1)\epsilon}}\sim\frac{1}{\sqrt{2N}}\;,$
(61)

where for the approximation we again used the natural value of $\epsilon$ and
$N\gg 1$. We observe that the coupling of $\tilde{a}_{H}$ remains roughly the
same as in the single axion case, while $\tilde{a}_{L}$ shows a suppression.
Overall, we can say that on top of the above mentioned mass splitting into two
physical states, their couplings differ by the factors $f_{H}$ and $f_{L}$.

For the axion-photon coupling, this modification results in

$\displaystyle g_{a\gamma}^{H/L}$ $\displaystyle=\frac{\alpha}{2\pi
f_{a}}\left(\frac{\mathcal{E}}{\mathcal{N}}-1.92\right)f_{H/L}(N,\epsilon)$
$\displaystyle=\frac{\alpha}{2\pi f_{a}}C_{a\gamma}f_{H/L}(N,\epsilon)\;,$
(62)

where $\mathcal{E}$ and $\mathcal{N}$ are the electromagnetic and QCD anomaly
coefficients, respectively, and $\alpha$ is the electromagnetic fine-structure
constant. In the second line, we defined the model dependent factor
$C_{a\gamma}$ for compactness. In order to bring this to a more useful form,
we express $f_{a}$ in terms of the physical mass states from (56), i.e.

$m_{a}=\frac{\sqrt{\Lambda_{\rm
QCD}^{3}m_{u}}}{f_{a}}=\begin{cases}\sqrt{1-\epsilon}\ m_{H}&:{\rm for}\
\tilde{a}_{H}\par\;,\\\ \sqrt{1+(N-1)\epsilon}\ m_{L}\;&:{\rm for}\
\tilde{a}_{L}.\end{cases}$ (63)

After plugging in (60) and (61), this results in

$\displaystyle g_{a\gamma}^{H}$
$\displaystyle=\frac{\alpha}{2\pi}\frac{C_{a\gamma}}{\sqrt{\Lambda_{\rm
QCD}^{3}m_{u}}}\sqrt{\frac{N-1}{N}}\ m_{H}$ (64) $\displaystyle
g_{a\gamma}^{L}$
$\displaystyle=\frac{\alpha}{2\pi}\frac{C_{a\gamma}}{\sqrt{\Lambda_{\rm
QCD}^{3}m_{u}}}\sqrt{\frac{1}{N}}\ m_{L}\;.$ (65)

For concreteness, let us investigate the axion-photon couplings in the case of
KSVZ axions. In Figure 5 we show the possible couplings for $N<10^{6}$. The
included projected sensitivities of future experiments is ranging into the
predicted band but most of the parameter space will not be covered in the near
future. Nevertheless, one has to keep in mind that the heavy axion, which
couples similarly to a scenario with just one axion (maximal enhancement is a
factor of $\sim 1.3$), has a realistic chance of being discovered in the near
future. In the framework of multiple axions, the next step after the discovery
of the ordinary axion would be to search for a second light weakly coupled
state. Even though one could think that the task is quite difficult because
the available parameter space is large, one has to keep in mind that there is
a clear functional dependence among the couplings and the masses between the
heavy and the light state. So imagine a situation that an axion has been
discovered and therefore $m_{H}$ and $g_{a\gamma}^{H}$ would be known. In the
expression for the coupling of the light axion (65), we can replace $m_{L}$ by
$m_{H}$ using (57). We can then eliminate $N$ by plugging in the expression
for the heavy axion (64), yielding an expression for $g_{a\gamma}^{L}$ in
terms of experimentally known parameters,

$\displaystyle g_{a\gamma}^{L}\sim\frac{\alpha
C_{a\gamma}m_{H}}{2\pi\sqrt{\Lambda_{\rm
QCD}^{3}m_{u}}}\sqrt{1-\left(\frac{2\pi g_{a\gamma}^{H}\sqrt{\Lambda_{\rm
QCD}^{3}m_{u}}}{\alpha C_{a\gamma}m_{H}}\right)^{2}}\;,$ (66)

This means that after measuring the coupling and the mass of the first axion,
the properties of the second are uniquely determined.

###  III.3 Modification of Axion Physics for one YM Sector

Another interesting case for kinetic mixing is the scenario of a single YM
sector with PQ scale $f_{a_{2}}$ and confinements scale $\Lambda_{\rm conf}$,
as discussed in Section II.5. The procedure of analyzing the kinetic mixing
remains unchanged and $N=2$ in this case.

The coupling to our photons after the step of choosing the proper kinetic
basis reads as

$g_{a\gamma}^{L/H}=\frac{\alpha}{2\pi
f_{a}}C_{a\gamma}\frac{1}{\sqrt{2}}\frac{1}{\sqrt{1\pm\epsilon}}\;.$ (67)

Since $\Lambda_{\rm conf}$ and $f_{a_{2}}$ differ from the analog quantities
in our sectors, the axion masses differ as well. This leads to a non-diagonal
mass matrix,

$\displaystyle\mathcal{L}_{\rm mass}$ $\displaystyle\sim
m_{a_{1}}^{2}a_{1}a_{1}+m_{a_{2}}^{2}a_{2}a_{2}$
$\displaystyle\sim\begin{pmatrix}\tilde{a}_{H}\\\
\tilde{a}_{L}\end{pmatrix}^{T}\begin{pmatrix}\frac{m_{a_{1}}^{2}+m_{a_{2}}^{2}}{1-\epsilon}&&\frac{m_{a_{1}}^{2}-m_{a_{2}}^{2}}{\sqrt{1-\epsilon}\sqrt{1+\epsilon}}\\\
\frac{m_{a_{1}}^{2}-m_{a_{2}}^{2}}{\sqrt{1-\epsilon}\sqrt{1+\epsilon}}&&\frac{m_{a_{1}}^{2}+m_{a_{2}}^{2}}{1+\epsilon}\end{pmatrix}\begin{pmatrix}\tilde{a}_{H}\\\
\tilde{a}_{L}\end{pmatrix}.$ (68)

As usual, this mass-term must be diagonalized by an orthogonal matrix with the
mixing angle

$\theta=\frac{1}{2}\arctan\left(\frac{m_{a_{1}}^{2}-m_{a_{2}}^{2}}{m_{a_{1}}^{2}+m_{a_{2}}^{2}}\frac{\sqrt{1-\epsilon^{2}}}{\epsilon}\right)\;,$
(69)

leading to the eigenvalues,

$\displaystyle
m_{H/L}=\frac{m_{a_{1}}^{2}+m_{a_{2}}^{2}\pm(m_{a_{1}}^{2}-m_{a_{2}}^{2})\sqrt{1+\frac{4\epsilon^{2}}{\left(\frac{m_{a_{1}}}{m_{a_{2}}}-\frac{m_{a_{2}}}{m_{a_{1}}}\right)^{2}}}}{1-\epsilon^{2}}\;.$
(70)

Expressing $a_{H}$ and $a_{L}$ in terms of the final eigenstates $A_{H}$ and
$A_{L}$, the axion-photon coupling from (67) becomes

$\displaystyle g_{A\gamma}^{H/L}$ $\displaystyle=\frac{\alpha}{4\pi
f_{a_{1}}}C_{a\gamma}\left(\frac{\cos{\theta}}{\sqrt{1+\epsilon}}\pm\frac{\sin{\theta}}{\sqrt{1-\epsilon}}\right)$
$\displaystyle\equiv\frac{\alpha}{4\pi f_{a_{1}}}C_{a\gamma}\kappa_{H/L}\;.$
(71)

Using this coupling and the masses, given in (70), the lifetime of the axions
is given by

$\displaystyle\tau(A_{H/L}\rightarrow
2\gamma)=\frac{2^{6}\pi}{\left(g_{A\gamma}^{H/L}\right)^{2}m_{H/L}}=\frac{2^{10}\pi^{3}}{\alpha^{2}C_{a\gamma}^{2}}\frac{f_{a_{1}}^{2}}{\kappa_{H/L}^{2}m_{H/L}^{3}}\;.$
(72)

For the axion to constitute a part of the dark matter,
$\tau(A_{H/L}\rightarrow 2\gamma)$ must be larger than the age of the
universe. This requirement translates to the viable region shown in Figure 6.
While we chose $\epsilon\sim 1/2$ in this figure, the viable region is
essentially independent of $\epsilon$.

Figure 6: Same as Figure. 3 but with axion kinetic mixing. The mixing results
in both axion states being able to decay into our photons. Requiring the
stability of the heavy axion state results in the viable region to lie below
the dashed orange line ($f_{a}\sim 10^{12}\,\mathrm{GeV}$) or below the dashed
purple line ($f_{a}\sim M_{\mathrm{P}}$).

The interpretation of this viable region is vastly unchanged from the
discussion in Section II.5. With kinetic mixing an additional upper bound of
the dark confinement scale emerges. Consequently, if the dark sector axion
makes up the dark matter, the confinement scale of the dark YM sector must be
in the range

$1\,\mathrm{eV}\lesssim\Lambda_{\rm conf}\lesssim 10^{12}\,\mathrm{GeV}\;.$
(73)

##  IV Conclusions

In this paper, we studied phenomenological consequences of adding axions to
dark YM sectors. These axions are required from quantum gravitational
arguments that promote the strong CP problem to a consistency problem. As
such, an axion becomes a necessary component not only in QCD but in every YM
group Dvali:2018dce ; dvali2022strong . We focused on two models, namely $N$
exact copies of the SM Dvali:2007hz ; Dvali:2007wp  and a single pure YM
sector.

We first required that the total axion density including the axions from all
sectors must not exceed the observed dark matter density. For $N$ exact SM
copies, using the misalignment mechanism in all sectors results in a limited
range of viable $N$ for a given PQ scale. Since the misalignment mechanism is
independent of the thermalization of the dark sectors, this applies even for
very large values of $N$ when the sectors are very dilute. For a single pure
YM sector, we find a relation between the dark confinement scale and the PQ
scale. In particular, requiring the dark sector axion to make up the dark
matter, the PQ scale does not exceed the Planck scale when the dark
confinement scale fulfills $\Lambda_{\rm conf}\gtrsim 1\,\mathrm{eV}$ in the
scenario with maximal misalignment. This minimal value is pushed higher for
smaller initial misalignment angles.

Furthermore, we showed that the contribution of additional axions to the
isocurvature fluctuations does not necessarily tighten the bound on
$H_{\mathrm{I}}$. In fact, both models exhibit regions in the parameter space
where the bound is essentially avoided. For $N$ exact SM copies the viable
parameter space requires either a moderate number of $N\sim 10^{2}$ copies or
a mechanism that allows for smaller PQ scales.

In models with several dark sectors, it is possible that particles from
different sectors collectively form structures with a mass spectrum suppressed
by $1/\sqrt{N}$ Dvali:2019ewm . While this phenomenon was originally studied
with fermions, we showed that it also applies to bosons with a repulsive self-
interaction in the regime of dilute boson stars. Since axions have attractive
self-interactions, this phenomenon cannot take place. Therefore, the presence
of multiple axions does not modify the mass spectrum of such axion stars.

Lastly, we studied non-gravitational communication between the dark sectors
from axion kinetic mixing. We showed that a mass splitting into $N-1$
degenerate states and one light state emerges. Interestingly, the latter is
lighter by a factor of $1/\sqrt{2}$ and has couplings weaker by a factor of
$1/\sqrt{N}$ with respect to the heavier axion states. For $N$ exact SM
copies, this allows us to completely predict the light axion once the heavier
one is discovered. For a single YM sector, the axion masses are different, so
that the axion from the dark sector could decay into photons. If this axion is
supposed to make up the dark matter, an upper bound on the dark confinement
scale appears by requiring the axion’s stability, i.e. $\Lambda_{\rm
conf}\lesssim 10^{12}\,\mathrm{GeV}$. This bound is for most parts of the
parameter space stronger than the bound $\Lambda_{\rm conf}\lesssim f_{a_{2}}$
arising from perturbative unitarity.

The present paper provides certain insights into the potential consequences of
multiple axions. As mentioned earlier, the situation with multiple axions
closely resembles the axiverse proposed by string theory Arvanitaki:2009fg .
However, in our case, the axions arise from quantum gravitational consistency
in theories with numerous hidden YM groups. Given that hidden YM groups and
multiple axions are predictions of string theory, exploring this connection
further is a valuable pursuit. This connection can in principle also be
applied to dark energy. If the dark energy is realized by quintessence
Ratra:1987rm ; Wetterich:1987fm , it would suggest the existence of a non-
trivial vacuum structure that needs to be eliminated by this scalar field. All
in all, it seems that axions could play a larger role in nature than initially
anticipated.

###### Acknowledgements.

We are grateful to Gia Dvali and Georg Raffelt for very useful discussions and
helpful comments on the manuscript. This work was partly supported by the
Sonderforschungsbereich SFB1258.

## References

  * (1) J. M. Pendlebury et al., “Revised experimental upper limit on the electric dipole moment of the neutron,” Phys. Rev. D 92 no. 9, (2015) 092003, arXiv:1509.04411 [hep-ex].
  * (2) J. R. Ellis and M. K. Gaillard, “Strong and Weak CP Violation,” Nucl. Phys. B 150 (1979) 141–162.
  * (3) G. Dvali, “Strong-$CP$ with and without gravity,” arXiv:2209.14219 [hep-ph].
  * (4) G. Dvali and C. Gomez, “Quantum Compositeness of Gravity: Black Holes, AdS and Inflation,” JCAP 01 (2014) 023, arXiv:1312.4795 [hep-th].
  * (5) G. Dvali, C. Gomez, and S. Zell, “Quantum Breaking Bound on de Sitter and Swampland,” Fortsch. Phys. 67 no. 1-2, (2019) 1800094, arXiv:1810.11002 [hep-th].
  * (6) G. Dvali, “$S$-Matrix and Anomaly of de Sitter,” Symmetry 13 no. 1, (2020) 3, arXiv:2012.02133 [hep-th].
  * (7) G. Dvali, C. Gomez, and S. Zell, “A Proof of the Axion?,” arXiv:1811.03079 [hep-th].
  * (8) C. Vafa and E. Witten, “Parity Conservation in QCD,” Phys. Rev. Lett. 53 (1984) 535.
  * (9) R. D. Peccei and H. R. Quinn, “Constraints Imposed by CP Conservation in the Presence of Instantons,” Phys. Rev. D 16 (1977) 1791–1797.
  * (10) R. D. Peccei and H. R. Quinn, “Constraints imposed by $\mathrm{CP}$ conservation in the presence of pseudoparticles,” Phys. Rev. D 16 (Sep, 1977) 1791–1797. https://link.aps.org/doi/10.1103/PhysRevD.16.1791.
  * (11) S. Weinberg, “A New Light Boson?,” Phys. Rev. Lett. 40 (1978) 223–226.
  * (12) F. Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40 (1978) 279–282.
  * (13) G. Dvali, “Three-form gauging of axion symmetries and gravity,” arXiv:hep-th/0507215.
  * (14) G. Dvali and G. R. Farrar, “Strong CP Problem with 10**32 Standard Model Copies,” Phys. Rev. Lett. 101 (2008) 011801, arXiv:0712.3170 [hep-th].
  * (15) G. Dvali, I. Sawicki, and A. Vikman, “Dark Matter via Many Copies of the Standard Model,” JCAP 0908 (2009) 009, arXiv:0903.0660 [hep-th].
  * (16) G. Dvali and M. Redi, “Phenomenology of 1032 Dark Sectors,” Phys. Rev. D80 (2009) 055001, arXiv:0905.1709 [hep-ph].
  * (17) G. Dvali, “Black Holes and Large N Species Solution to the Hierarchy Problem,” Fortsch. Phys. 58 (2010) 528–536, arXiv:0706.2050 [hep-th].
  * (18) G. Dvali and M. Redi, “Black Hole Bound on the Number of Species and Quantum Gravity at LHC,” Phys. Rev. D77 (2008) 045027, arXiv:0710.4344 [hep-th].
  * (19) S. I. Blinnikov and M. Y. Khlopov, “ON POSSIBLE EFFECTS OF ’MIRROR’ PARTICLES,” Sov. J. Nucl. Phys. 36 (1982) 472.
  * (20) E. W. Kolb, D. Seckel, and M. S. Turner, “The Shadow World,” Nature 314 (1985) 415–419.
  * (21) D. J. Gross, J. A. Harvey, E. J. Martinec, and R. Rohm, “The Heterotic String,” Phys. Rev. Lett. 54 (1985) 502–505.
  * (22) L. J. Dixon, J. A. Harvey, C. Vafa, and E. Witten, “Strings on Orbifolds,” Nucl. Phys. B 261 (1985) 678–686.
  * (23) O. Lebedev, H. P. Nilles, S. Raby, S. Ramos-Sanchez, M. Ratz, P. K. S. Vaudrevange, and A. Wingerter, “A Mini-landscape of exact MSSM spectra in heterotic orbifolds,” Phys. Lett. B 645 (2007) 88–94, arXiv:hep-th/0611095.
  * (24) V. Braun, Y.-H. He, B. A. Ovrut, and T. Pantev, “A Heterotic standard model,” Phys. Lett. B 618 (2005) 252–258, arXiv:hep-th/0501070.
  * (25) G. Dvali, “Entropy Bound and Unitarity of Scattering Amplitudes,” JHEP 03 (2021) 126, arXiv:2003.05546 [hep-th].
  * (26) J. Halverson, C. Long, A. Maiti, B. Nelson, and G. Salinas, “Gravitational waves from dark Yang-Mills sectors,” JHEP 05 (2021) 154, arXiv:2012.04071 [hep-ph].
  * (27) W.-C. Huang, M. Reichert, F. Sannino, and Z.-W. Wang, “Testing the dark SU(N) Yang-Mills theory confined landscape: From the lattice to gravitational waves,” Phys. Rev. D 104 no. 3, (2021) 035005, arXiv:2012.11614 [hep-ph].
  * (28) G. Dvali, “Nature of Microscopic Black Holes and Gravity in Theories with Particle Species,” Int. J. Mod. Phys. A 25 (2010) 602–615, arXiv:0806.3801 [hep-th].
  * (29) M. Ettengruber, “Neutrino physics in TeV scale gravity theories,” Phys. Rev. D 106 no. 5, (2022) 055028, arXiv:2206.00034 [hep-ph].
  * (30) G. Dvali and O. Pujolas, “Micro Black Holes and the Democratic Transition,” Phys. Rev. D 79 (2009) 064032, arXiv:0812.3442 [hep-th].
  * (31) G. Dvali, E. Koutsangelas, and F. Kuhnel, “Compact Dark Matter Objects via $N$ Dark Sectors,” Phys. Rev. D 101 (2020) 083533, arXiv:1911.13281 [astro-ph.CO].
  * (32) J. Preskill, M. B. Wise, and F. Wilczek, “Cosmology of the Invisible Axion,” Phys. Lett. B 120 (1983) 127–132.
  * (33) M. Dine and W. Fischler, “The Not So Harmless Axion,” Phys. Lett. B 120 (1983) 137–141.
  * (34) L. F. Abbott and P. Sikivie, “A Cosmological Bound on the Invisible Axion,” Phys. Lett. B 120 (1983) 133–136.
  * (35) G. R. Dvali, “Removing the cosmological bound on the axion scale,” arXiv:hep-ph/9505253.
  * (36) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, “String Axiverse,” Phys. Rev. D 81 (2010) 123530, arXiv:0905.4720 [hep-th].
  * (37) K. J. Mack and P. J. Steinhardt, “Cosmological Problems with Multiple Axion-like Fields,” JCAP 05 (2011) 001, arXiv:0911.0418 [astro-ph.CO].
  * (38) G. Dvali, “Bounds on quantum information storage and retrieval,” Phil. Trans. A. Math. Phys. Eng. Sci. 380 no. 2216, (2021) 20210071, arXiv:2107.10616 [hep-th].
  * (39) G. Dvali, “Saturon Dark Matter,” arXiv:2302.08353 [hep-ph].
  * (40) A. R. Zhitnitsky, “On Possible Suppression of the Axion Hadron Interactions. (In Russian),” Sov. J. Nucl. Phys. 31 (1980) 260.
  * (41) M. Dine, W. Fischler, and M. Srednicki, “A Simple Solution to the Strong CP Problem with a Harmless Axion,” Phys. Lett. B 104 (1981) 199–202.
  * (42) J. E. Kim, “Weak Interaction Singlet and Strong CP Invariance,” Phys. Rev. Lett. 43 (1979) 103.
  * (43) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, “Can Confinement Ensure Natural CP Invariance of Strong Interactions?,” Nucl. Phys. B 166 (1980) 493–506.
  * (44) L. Di Luzio, F. Mescia, and E. Nardi, “Window for preferred axion models,” Phys. Rev. D 96 no. 7, (2017) 075003, arXiv:1705.05370 [hep-ph].
  * (45) V. Plakkot and S. Hoof, “Anomaly ratio distributions of hadronic axion models with multiple heavy quarks,” Phys. Rev. D 104 no. 7, (2021) 075017, arXiv:2107.12378 [hep-ph].
  * (46) J. Diehl and E. Koutsangelas, “Dine-Fischler-Srednicki-Zhitnitsky-type axions and where to find them,” Phys. Rev. D 107 no. 9, (2023) 095020, arXiv:2302.04667 [hep-ph].
  * (47) O. Sakhelashvili, “Consistency of the dual formulation of axion solutions to the strong CP problem,” Phys. Rev. D 105 no. 8, (2022) 085020, arXiv:2110.03386 [hep-th].
  * (48) R. Foot, H. Lew, and R. R. Volkas, “A Model with fundamental improper space-time symmetries,” Phys. Lett. B 272 (1991) 67–70.
  * (49) C. G. Callan, Jr., R. F. Dashen, and D. J. Gross, “Toward a Theory of the Strong Interactions,” Phys. Rev. D 17 (1978) 2717.
  * (50) D. J. Gross, R. D. Pisarski, and L. G. Yaffe, “QCD and Instantons at Finite Temperature,” Rev. Mod. Phys. 53 (1981) 43.
  * (51) P. Fox, A. Pierce, and S. Thomas, “Probing a QCD String Axion with Precision Cosmological Measurements,” arXiv e-prints (Sept., 2004) hep–th/0409059, arXiv:hep-th/0409059 [astro-ph].
  * (52) Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * (53) G. G. Raffelt, “Astrophysical axion bounds,” Lect. Notes Phys. 741 (2008) 51–71, arXiv:hep-ph/0611350.
  * (54) F. Takahashi, W. Yin, and A. H. Guth, “QCD axion window and low-scale inflation,” Phys. Rev. D 98 no. 1, (2018) 015042, arXiv:1805.08763 [hep-ph].
  * (55) E. Koutsangelas, “Removing the cosmological bound on the axion scale in the Kim-Shifman-Vainshtein-Zakharov and Dine-Fischler-Srednicki-Zhitnitsky models,” Phys. Rev. D 107 no. 9, (2023) 095009, arXiv:2212.07822 [hep-ph].
  * (56) D. E. Kaplan and R. Rattazzi, “Large field excursions and approximate discrete symmetries from a clockwork axion,” Phys. Rev. D 93 no. 8, (2016) 085007, arXiv:1511.01827 [hep-ph].
  * (57) T. Kobayashi, R. Kurematsu, and F. Takahashi, “Isocurvature Constraints and Anharmonic Effects on QCD Axion Dark Matter,” JCAP 09 (2013) 032, arXiv:1304.0922 [hep-ph].
  * (58) I. I. Tkachev, “On the possibility of Bose star formation,” Phys. Lett. B 261 (1991) 289–293.
  * (59) H. Zhang, “Axion Stars,” Symmetry 12 no. 1, (2019) 25, arXiv:1810.11473 [hep-ph].
  * (60) P.-H. Chavanis, “Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions: I. Analytical results,” Phys. Rev. D 84 (2011) 043531, arXiv:1103.2050 [astro-ph.CO].
  * (61) L. Berezhiani, G. Cintia, and M. Warkentin, “Core fragmentation in simplest superfluid dark matter scenario,” Phys. Lett. B 819 (2021) 136422, arXiv:2101.08117 [astro-ph.CO].
  * (62) L. Berezhiani, G. Cintia, and J. Khoury, “Thermalization, fragmentation, and tidal disruption: The complex galactic dynamics of dark matter superfluidity,” Phys. Rev. D 107 no. 12, (2023) 123010, arXiv:2212.10577 [astro-ph.CO].
  * (63) R. Foot, “Mirror dark matter: Cosmology, galaxy structure and direct detection,” Int. J. Mod. Phys. A 29 (2014) 1430013, arXiv:1401.3965 [astro-ph.CO].
  * (64) C. O’Hare, “cajohare/axionlimits: Axionlimits.” https://cajohare.github.io/AxionLimits/, July, 2020.
  * (65) G. Dvali, E. Koutsangelas, and F. Kuhnel, “Compact Dark Matter Objects via $N$ Dark Sectors,” Phys. Rev. D 101 (2020) 083533, arXiv:1911.13281 [astro-ph.CO].
  * (66) B. Ratra and P. J. E. Peebles, “Cosmological Consequences of a Rolling Homogeneous Scalar Field,” Phys. Rev. D 37 (1988) 3406.
  * (67) C. Wetterich, “Cosmology and the Fate of Dilatation Symmetry,” Nucl. Phys. B 302 (1988) 668–696, arXiv:1711.03844 [hep-th].