Exploring the Universe with Dark Light Scalars

We study the cosmology of a dark sector consisting of (ultra) light scalars. Since the scalar mass is radiatively unstable, a special explanation is required to make it much smaller than the UV scale. There are two well-known mechanisms for the origin of scalar mass. The scalar can be identified as a pseudo-Goldstone boson, whose shift symmetry is explicitly broken by non-perturbative corrections, like the axion. Alternatively, it can be identified as a composite particle like the glueball, whose mass is limited by the confinement scale of the theory. In both cases, the scalar can be naturally light, but interaction behavior is quite different. The lighter the axion (glueball) mass is, the weaker (stronger) the interaction strength is. We consider the dark axion whose shift symmetry is anomalously broken by the hidden non-abelian gauge symmetry. After the confinement of the gauge group, the dark axion and the dark glueball get masses and both form multicomponent dark matter. We carefully consider the effects of energy flow from the dark gluons to the dark axions and derive the full equations of motion for the background and the perturbed variables. The effect of the dark axion-dark gluon coupling on the evolution of the entropy and the isocurvature perturbations is also clarified. Finally, we discuss the gravo-thermal collapse of the glueball subcomponent dark matter after the halos form, in order to explore the potential to contribute to the formation of seeds for the supermassive black holes observed at high redshifts. With simplified assumptions, the glueball subcomponent dark matter with the mass of $\sim 0.01-1{\rm MeV}$, and the axion dark matter with the GUT scale decay constant (the mass of ${\cal O}(10^{-18})\,{\rm eV}$) can provide the hint on the origin of the supermassive black holes at high redshifts.


I. INTRODUCTION
The discovery of the Higgs boson completed the Standard Model which explains most of the phenomena in the Universe including nuclei and atoms and their interactions. However, various studies based on precise measurements of the cosmic microwave background (CMB), velocity distributions of stars and galaxies, and largescale structures (LSS) show that only 4% of the Universe is understood by our knowledge of the Standard Models and 96% must be filled with dark matter and dark energy that we do not know their origin [1].
Among dark matter candidates, particles classified as 'weakly interacting massive particle' (WIMP) have been considered as the best candidate. This is because their freeze-out relic abundance can naturally explain the present amount of dark matter, can be tested by different searches, and their mass scale can be related to new physics that explains why the electroweak scale is stable against various quantum corrections. However, there is no conclusive evidence of WIMP dark matter so far, and the alternative candidates got more attention in recent years. The masses of these candidates are not limited to * whqnrjs@snu.ac.kr † htkim428@snu.ac.kr ‡ hdkim@phya.snu.ac.kr § csshin@ibs.re.kr a narrow range, and various ideas have been proposed particularly focusing on the detection possibility [2]. For such a broad range of dark matter mass, no clear guiding principle exists to specify a natural range of it. Especially, it is extremely unnatural to consider a light scalar unless there is a special reason for it. Although theoretical explanations are required, it is quite interesting to note that an ultra-light scalar dark matter is allowed by cosmological/astrophysical observations, because the scalar can act as an oscillating classical field whose averaged equation of state is the same as that of the cold dark matter (CDM). Interesting phenomena also can arise from the field nature of the scalar. The fuzzy dark matter [3] and the QCD axion [4][5][6][7] are good examples. Both candidates, the axion or axion like particle, could be naturally light by the approximate shift symmetry of the axion, which is explicitly broken by the controllable non-perturbative corrections. The consequence of the global symmetry is that the lighter the axion is, the weaker its interaction is. Therefore, the large occupation number of the (ultra) light axion is allowed, and it can be described by the evolution of scalar condensate.
On one hand, there is another way to obtain a light scalar. When the asymptotically free gauge group has confinement at low energy, the gauge fields are confined into the scalar particles, the glueball, whose mass is limited by the confining scale. Unlike in the case of the axion, the lighter the glueballs is, the stronger the scattering cross-section among glueballs is. The number chang-arXiv:2010.10880v1 [hep-ph] 21 Oct 2020 ing interactions are active, so the occupation number is always limited by its temperature. In this case, the relic density is determined by the freeze-out mechanism. Because of this property, the light dark glueball becomes a good candidate for self-interacting dark matter (SIDM) [8,9], which is one of the ways to make the cored density profile around the center of galaxies [10]. As a subcomponent dark matter, it can also play the role to suppress small scale perturbation [11]. These two mechanisms to obtain a light scalar dark matter provide completely different microscopic nature of dark matter, and yield different predictions for small scale structure.
In this paper we study the cosmology of light scalar dark matter, focusing on the origin of its mass and the consequences associated with it. We cover two mechanisms discussed above at the same time in a minimal set-up: the dark sector consisting of the axion and the confining hidden gauge symmetry without light fermions. In this set-up, the axion scalar potential is dynamically generated by a Chern-Simons type coupling between the axion and the dark gauge field, gluon. As the gauge symmetry confines, the axion and the glueball get masses and become a part of multicomponent dark matter. The idea that the mass of the ultra-light scalar dark matter originates from a confining hidden gauge symmetry was studied in [12,13], but there is no comprehensive study for its cosmological evolution.
We derive a complete set of the equations of motion for the background and the perturbed variables of the dark axion and the dark gluon/glueball densities. Since the axion-gluon coupling provides energy transfer from the dark gluon to the dark axion, we clarify its effect on thermodynamics of the dark gluon fluid and entropy evolution. We also quantify the transfer of the isocurvature perturbation through the same coupling from each component to the other.
The multicomponent dark matter which simultaneously contains feebly interacting particles and strongly interacting particles has an interesting cosmological consequence. If the glueball becomes a subcomponent dark matter, the strong self-interaction between the glueballs can help form the seed of the supermassive black holes [14,15]. This can provide a possible answer to the question on the origin of observed supermassive black holes at high redshifts [16][17][18][19]. We discuss the parameter space to provide the solution and possible caveats.
The paper is organized as follows. In Sec. II, we establish basic formalism from the Lagrangian to the dynamical equations of the background and perturbation variables of the coupled axion-gluon fluid. In Sec. III, we set up the model for the evolution of the glueballs and axion as dark matter. The parametric dependence of the relic abundance of the glueballs and the axion is also presented. Sec. IV is devoted to the evolution of the perturbed variables. For the initial conditions, we have three modes: adiabatic perturbation, isocurvature perturbation induced by the initial misalingnment of the axion and by the temperature fluctuation of the dark gluon fluid. In Sec. V, the implications of the glueball subcomponent dark matter for the early formation of the supermassive black holes are discussed. Sec. VI is discussion.

II. AXION DARK MATTER AND CONFINING DARK SECTOR
A. General description of the model Our starting Lagrangian of dark sector is composed of the ultra-light axion φ whose field range is 2πf a , and the dark gauge symmetry with the confinement scale Λ (Λ f a ). The coupling between the axion and the dark gauge bosons are given as where G a µν is the dark gluon field strength and g denotes the dark gauge coupling. For illustration, we take SU (N ) as the dark gauge group in our set-up. Below the confinement scale, the dynamics of the gauge fields can be described by the composite bosons, the glueballs. Considering the lightest glueball ϕ g , the effective Lagrangian of ϕ g can be expanded in (4π/N )(ϕ g /m g ) [20][21][22][23] where m g is the glueball mass of O(Λ). The axion also gets a scalar potential by the gluodynamics, which can be written as a power series in (φ/N f a ) 2 around its minimum (φ = 0) [24,25] In this expansion, the axion mass is given by where V 4 is the four-dimensional volume. The 1/N 2 factor for the quartic term of the axion potential leads to the suppression of the anharmonic effects as long as the initial misalignment of the axion field is φ i f a . The glueball and the axion can be naturally light by lowering the scale of confinement. The difference between the two is that as their masses get smaller, the glueballs interact more strongly with each other, whereas the axions interact more feebly due to the additional suppression factor Λ/f a 1. The glueball is not a lightest particle in our set-up, so it is unstable. In the lattice calculation [26,27], the glueball mass depends on the axion field value as (2.5) Through the quadratic and cubic interactions of the glueball, we get the leading interaction for the decay of the glueball to two axions, From Eq. (2.6), the life-time of the glueball is estimated as τ ϕg ∼ 10 17 Gyr (f a /10 13 GeV) 4 (GeV/m g ) 5 . In the parameter space we will focus on, the glueball is cosmologically stable, so that both axion and glueball are dark matter of the Universe. For cosmology, we consider the case that the dark sector and the visible sector are thermally decoupled from the beginning. Therefore, dark gluons/glueballs are thermalized by their own interactions at a temperature T g , which is different from the photon temperature T γ . Starting from the gluon fluid at high temperatures, as T g drops and crosses the dark critical temperature T g,c = O(Λ), the confining phase transition occurs. Below T g,c , all gluons are confined into the glueballs, and the evolution is described by the massive glueball fluid. The dark gluon temperature also affects the evolution of the dark axions. The leading term of the axion potential induced by the gluo-thermodynamics is where the axion mass m a (T g ) is well described by the dilute instanton gas approximation in the deconfining phase, with η a = 11N/6 − 2. For N = 3(4), η a = 3.5(5.3) [28]. After the confinement (T g T g,c ), the axion mass is saturated as its zero temperature value, m a (T g ) m a [29].
Actually, the temperature dependence of the axion potential implies the existence of the energy flow from the gluon fluid to the axions as the temperature decreases. Then, a natural question is whether or not the entropy of the dark gluon also evolves during the energy transfer. In order to make it clearer, let us discuss the gluothermodynamics in more detail. The free energy density of the gluon/glueball fluid can be evaluated from the gluon partition function for a given temperature T g . If the topological θ-term (θ ≡ φ/f a ) is vanishing, the free energy density f g is only the function of the temperature as f g (T g ) = −p g , where p g is the pressure of the gluon/glueball fluid. The energy (ρ g ) and entropy (s g ) densities are obtained by the thermodynamic relations, s g = −df g /dT g and ρ g = T g s g − p g . On the other hand, the situation is a little bit different for the non-vanishing θ-term. Because the gluon partition function also depends on θ, the free energy density (the negative of the pressure for the gluon fluid with the θ-term) is evaluated as [26] f g+θ (T g ) = −p g + V (T g , φ), (2.9) where p g ≡ −f g (T g ). The second term of the RHS represents the contribution of the vacuum energy density generated by the non-perturbative gluo-thermodynamics. The energy density of the gluon fluid with the θ-term also can be decomposed into the sum of the vacuum energy density and the pure gluonic contribution as ρ g+θ (T g ) = ρ g + V (T g , φ). Then, the thermodynamics relations provide (2.10) As a result, the entropy and the energy density of the gluon/glueball fluid depend not only on the temperature but also on the axion field value as Eq. (2.10) and From the continuity equation of the dark sector, we can explicitly show that the entropy of the dark sector in a comoving volume s g a 3 is conserved for whatever value of φ and its evolution except the period of the confining phase transition. Before discussing the explicit evolution of each component, for completeness, we address general equations of motion for axion and gluon as the fluids including their homogenous and perturbation parts.

B. Dynamics of the axion-gluon/glueball fluids
The axion dark matter is described by the evolution of the classical fields, φ(x). The dark gluon/glueball densities and their perturbations can be parameterized by its temperature evolution T g (x) and φ(x) as discussed in the previous section. In this context, {φ(x), T g (x)} are good variables to derive full equations of motion of dark sector including their perturbations. Considering the fluid description, the evolution of energy densities and pressures are deduced from the evolution of φ and T g with the help of the Einstein equations, gluo-thermodynamics and the lattice calculation.
In order to include perturbations, we introduce the gauge, the conformal Newtonian gauge, for the inhomogeneous part of the metric tensor for the time derivatives. In Sec. II and IV, the conformal time is mostly taken as the argument of the time dependent variables, whereas the proper time t is used for the discussion about the late time effects of dark matter in Sec. V. Expanding φ and T g near the homogeneous solutions, the equations of motion of the background axion field is given by The corresponding energy density and pressure are w a ≡ p a /ρ a is the equation of state for the axion. Since the axion-gluon fluid is isolated from the visible sector, the total energy and pressure of dark sector should satisfy the continuity equation without source terms ρ a+g + 3H(ρ a+g + p a+g ) = 0. (2.19) This leads to the evolution of the gluon/glueball fluid as where w g = p g /ρ g . Together with Eq. (2.10), (2.11), we can derive the conservation of the dark entropy as we claimed, For the evolution of the perturbed variables, the Fourier mode of the axion field δφ with the wavenumber k is evolved as The axion fluid perturbation variables δρ a , δp a , v a and π a are derived as The definition and convention of the variables in above equations follow Ref. [30][31][32]. From Eq. (2.11), the gluon/glueball fluid perturbations δp g and δρ g are also related with δT g and δφ as Similarly, from Eqs. (2.22) and (2.23), we derive the equations of motion for the fluid perturbations.
where δ a = δρ a /ρ a , u a = (1+w a )v a , and same definitions for δ g , u g . Note that {δp a , δp g , δT g , δφ} can be expressed by {δ a , δ g , u a } using Eq. (2.23) and (2.24), and solve the equations of {δ a , δ g , u a , u g }. Otherwise, the full equations also can be expressed by {δφ, δT g , u g }. Compared with the isolated axion and gluon/glueball perturbations, the equations contain the terms proportional to ∂V /∂T g , which originate from the non-perturbative interactions between the gluon and the axion. With the T g dependent axion potential given by Eq. (2.8), the contribution of these terms becomes larger as T g approaches to T g,c .
The amount of dark gluons can be parameterized by the ratio between the entropies of the dark sector and the visible sector, s g /s SM . Although dealing with the entropy ratio between two sectors would be easier to trace the evolution of the densities, in order to have a more intuitive picture about how cold the gluons are compared to the SM sector, we will use the ratio parameter between the temperatures. Taking a period around the phase transition, the photon temperature when T g arrives at T g,c is denoted by T γ,c . Then, we define the ratio parameter r as where s g (s SM ) is the entropy density of the gluon fluid (the SM sector), and g * S is the effective number of degrees of freedom in entropy for the SM sector.
Up to now, we have ignored the effect of dissipation induced by the background dark gluon plasma. Although, it is not crucial in our discussion, let us clarify the condition that our assumption is invalid. Including the friction term (γ fr ) for the axion's motion, the equation of motion of the background axion is written as [33] φ In the deconfining phase, γ fr = Γ sph (T g )/f 2 a T g , where the sphaleron rate is estimated as [34] for g 2 /4π 0.1. Around the critical temperature, T g ∼ T g,c the gauge coupling can be as large as g 2 N/4π = O(1). In this regime, the sphaleron rate is expected as Γ sph (T g,c ) ∼ T 4 g,c from the dimensional analysis. After the confining phase transition, no reliable calculation has been done so far. A crude estimation based on the dimensional analysis is that the dissipation rate is at most proportional to the entropy density (or number density) of the glueballs as γ fr ∼ s g /f 2 a . Because the time dependence of the Hubble rate and the dissipation rate are given as H ∝ a −2 (a −3/2 ) in radiation dominated era (matter dominated era) and γ fr ∝ a −3 , the gluon induced friction term is important only when the temperature of the visible sector becomes greater than Comparing this with the temperature when the axion starts to oscillate (T γ < TeV), it is always irrelevant.

A. Evolution of the background gluon and glueballs
In Eqs. (2.18-2.25), we establish the continuity equations for the background and perturbative variables based on the equations of motion of the axion field, and gluo-thermodynamics. On one hand, the dependence on the temperature of ρ g and p g should be clarified. It is particularly easy for T g T g,c and T g T g,c . In the limit of T g T g,c , gluons are just relativistic fluid. Therefore, and the axion-gluon interaction is negligible because However, as the gluon temperature approaches T g,c , Eq. (3.1) does not hold. The true dependence on the gluon temperature can only be figured out by the lattice calculation [35][36][37]. We adopt the lattice data with θ = 0, and deduce the densities for nonzero θ case based on the data of gluon/glueball pressure p g (T g ). At T g,c , the dark gluons are combined into dark glueballs with masses of multiples of the confining scale Λ [38][39][40][41]. Here, Λ is defined with a certain regularization scheme. Taking the M S scheme, the lattice results with the functional method [42,43] show that The phase transition is first-order if N > 2. In order to understand how long the phase transition happens, we compare the energy density in the deconfined phase with that of the confined phase at T g,c . The former is the energy density of the gluon fluid of O(0.1N 2 T 4 g,c ), and the latter is the sum of the tower of the glueballs of O(0.1T 4 g,c ). Thus, as N increases, larger latent heat is released and the transition period becomes longer. Let us shortly discuss how we evaluate the energy density of the glueballs.
Since the lightest glueball mass is calculated as [39,44] all glueballs are non-relativistic at T g,c . With the help of the spectral densityρ(m), the energy density of glueballs can be written as ρ(m) is approximated by the sum of the discrete low-lying resonances with masses m g (J P C ), where J (angular momentum), P (parity), and C (charge conjuration) are the eigenvalues of the states, and the continuum spectrum of the Hagedorn tower where n 2 = 1, n N ≥3 = 2, and the threshold mass m th is usually chosen as 2m g . The Hagedorn temperature T H is related with T g,c as [45] T Note that even if N increases, the number of degrees of freedom of the glueballs does not increase. Therefore, the contribution from the Hagedorn glueballs is insensitive to N and becomes O(0.01 − 0.1)T 4 g,c . The low-lying glueball contributions are O(10 − 50)% of it. We adopt the spectrum of the glueballs in Ref. [39].
The lower limit of the actual transition period is given by the period assuming the quasi-equilibrium transition, i.e. the pressures of deconfining/confining phases are equal and the latent heat is released adiabatically as the Universe expands. In this case, the temperatures of the confining and deconfining phases are the same and maintained at around T g,c . The duration of the phase transition is estimated by the conservation of the dark entropy: where a cf (a ci ) is the scale factor when the phase transition ends (starts), s gluon (T g,c ) (s glueball (T g,c )) denotes the entropy density of the dark gluon (glueball) at T g = T g,c .
The N -dependence of the duration is obtained from . The phase transition becomes stronger first-order as N increases. Thus, the additional entropy is generated during the transition and makes the glueballs hotter than the previous estimation. However, this effect is negligible unless N 4π, because the nucleation temperature is just around T g,c and the strong interactions of the gluon and glueball fluids provide a large friction coefficient for the bubble wall propagation. For N = 3, a cf 2a ci is obtained numerically, which is well matched with our parametric estimation N 2/3 . The assumption of dark entropy conservation will be kept in the following discussion.
After the transition, the Hagedorn contribution becomes negligible as the scale factor increases by an order of magnitude compared to a cf . The contribution of the low-lying glueballs maintains chemical equilibrium by two-to-two and three-to-two scatterings whose rates are estimated as where v f is the relative velocity of the final particles from the scattering. The detailed description of the evolution of low-lying stable glueballs is provided in [46]. If the P and T violating topological θ-term is absent, these glueballs can be classified as the eigenstates of J P C . In our set-up, because of the axion's motion, the time dependent mixings can arise between different P eigenstates, e.g. between 0 ++ and 0 −+ . As a result, the terms like may be generated. It turns out that theses mixings do not lead to the decay of heavier state to the light states by the kinematic reason. We expect that the stability of low lying glueballs is not altered by the axion's evolution.
As the universe expands, the 3 → 2 process freezes out when the most of the 2 → 2 processes still active. This is because the interaction rate of the 3 → 2 process is proportional to the square of the number density of the glueballs, while that of the 2 → 2 processes is linearly proportional to the number density of the glueballs. Therefore, the relative chemical equilibrium holds between different glueball states. Since the 3 → 2 process of the lightest glueball keeps the chemical potential of the lightest glueball to be zero, the chemical potential of all glueball states is zero until the 3 → 2 process of the lightest glueball freezes out. Before the freeze-out of 3 → 2 interactions, a single glueball state with a mass m g in thermal equilibrium is given by (3.11) The total energy density and pressure of the glueballs can be obtained by summing Eq. (3.11) over all lowlying stable glueballs. It is found that the contributions other than the lightest glueball become negligible for T g 0.5T g,c .
A distinguishing properties of the dark matter with the three-to-two self-interaction is the scaling behavior as the Universe expands. The temperature of the glueballs drops much slower than that of the photon during its chemical equilibrium maintaining s g ∝ 1/a 3 , so T g ∼ 1/ ln a. The energy density drops faster than that of a cold dark matter, since the three-to-two self-interaction converts the mass energy to the kinetic energy. This behavior ends when the process freezes out at T g = T g,f o with After that, the glueballs act as free-streaming particles. Using Eq. (3.13) and the dark entropy conservation, the An example of the evolution of the energy density for the case of the axion-dominated (left) and the glueballdominated (right). After the confinement, the number-changing self-interaction of the glueballs reduces its total number only logarithmically. In the axion-dominated scenario, the energy density of the axion is O ρg/N 2 at the time of confinement and dominates that of the glueballs afterward. In the glueball-dominated scenario, the energy density of the gluon and the glueballs is much larger than that of the axion for all epoch. Tg,c/r is the photon temperature when the confining phase transition of the dark sector starts. For the region above the line R(r, fa) = 1, the oscillation of the axion starts earlier than the confinement phase transition and the gluebll is the dominant dark matter component with the mass mg 6 rTγ,c, while below R(r, fa) = 1, the axion starts to oscillate after the transition. Here, we did not impose the constraints from the current bound, which are discussed in text.
freeze-out temperature is evaluated as if it happens during radiation dominated era. As a specific example, for N = 3, r = 0.01, and m g = 1 MeV, we get T g,f o 0.04 m g 0.2 T g,c . (3.15) The relation T g,f o = O(0.2)T g,c is not much sensitive to the values of r and m g that we are interested in. During the evolution of the glueballs, the photon temperature also evolves. When the dark glueballs freezeout, the photon temperature becomes We can also easily evaluate the case that the freeze-out of the dark glueball happens after the matter-radiation equality for m g < keV.
So far, we have specified all history of the gluons and glueballs in order to identify the time dependence of the glueball temperature T g (τ ) for Eq.  The lattice studies provide the coefficient of each term in perturbation expansion of the axion scalar potential Eq. (2.3) as [47].
One interesting feature of the axion potential is that it is not a single branch, but multiple (N ) branches, where for each branch the period of the scalar potential is 2πN f a [48]. The general expression of the scalar potential for a kth branch is where h(ψ) is the 2π-periodic function. At present, calculation of the full shape of the potential is not available.
Besides the lattice study, the scalar potential for the part of axion field range was studied in the large N limit using the holographic description of the pure SU (N ) gauge theory [49,50]. The 't Hooft coupling λ = g 2 N at the KK scale in the dual gravity theory can be matched with that deduced from the axion potentials of the lattice calculation Eq. (3.18). We find that λ = 10 − 20 gives a reasonable matching. At high temperatures of the gluons, the instanton approximation for the axion potential is valid, and there is a single branch. During the phase transition, branches will emerge, and the axion can be located in a different branch in a different patch of the Universe. If each branch provides a stable axion trajectory, we have to consider the effect of them seriously.
Following the approach of the holographic description as in [49], we estimate the tunneling rate between kth to k − 1th branches.
where the Euclidean action is This can be significantly large only when N 10 3 . Therefore, in our consideration, all branches with higher energy densities are quite unstable, and the transition to the lowest energy state will occur almost immediately. As a result, the effective potential of the axion is well described by (3.22) and one can think the evolution of the axion within the range 2πf a . Without worrying about the effect of other branches, Eq. (2.16) gives for φ f a . If the second term of the LHS is much larger than the third term, the axion field is approximately constant because of the large Hubble friction. This is the slow-roll limit. In the opposite case, the axion field oscillates with the oscillation frequency m a (T g ). Such evolution can be well approximated by the simple transition at a = a osc to satisfy 3H = am a (T g ) (3H = m a (T g )).
For each epoch, where θ i is the initial misalignment angle of the axion field, A(τ ) is slowly varying function with A /A am a (T g ). The axion acts like dark energy during a < a osc , while for a > a osc , the axion plays the role of cold dark matter because w a 0 by averaging out the fast oscillation. The initial axion value φ i is not deterministic. Since both θ i 1, and |θ i − π| 1 need some tuning or special model building, here we take Since the axion's mass depends on the history of the dark gluons (Eq. (2.8)), there are two characteristic scales which determine the evolution history of the axion: a osc (onset of the axion oscillation) and a ci (onset of the confining phase transition). As the scale factor approaches a ci , the contribution of the gluons to the axion's potential becomes substantial, and the axion mass is saturated. The evolution of the axion mass is smooth unless N is very large. There is no significant distinction between the deconfined and confined phases.
With the help of the lattice data, we find that the following quantity determines whether or not the axion starts to oscillate before the confining transition. If R(r, f a ) > 1, the axion starts to oscillate before the phase transition. The corresponding photon temperature becomes For the evolution of the axion density, the initial density of the axion at T γ = T γ,osc is approximated as After that, the axion field will oscillate with the time dependent frequency. For the combination Eq. (2.18) gives Since the fast oscillation of the axion field means w a = 0, N a 0 and N a is nearly conserved, i.e. ρ a /m a (T g ) ∝ 1/a 3 . Therefore, if R(r, f a ) > 1, the axion starts to oscillate before the confining phase transition (T γ,osc > T γ,c ), and the present relic density of the axion dark matter becomes (3.32) If R(r, f a ) < 1, the axion oscillates after the confining phase transition (T γ,osc < T γ,c ). The corresponding axion dark matter density is given by (3.33) As shown in Eq. (3.32), the glueballs dominate the total energy density of dark matter if the axion oscillation happens earlier than the confining phase transition. The reason is simply that initial axion energy density before the oscillation is bounded by the confining scale Λ 4 . On on hand, if r 2 /f a becomes small enough, so that the axion starts to oscillate after the phase transition, the axion becomes a dominant component of dark matter. Fig. 1 shows the evolution of dark matter densities for both glueball and axion domination scenarios. Fig. 2 shows the parametric dependence of the current relic abundance of the dark matter. The axes are represented by T γ,c , the photon temperature when the confining phase transition of the dark gauge sector starts, and m a , the zero temperature axion mass, defined in Eq. (2.4) with Eq. (3.18), respectively. As we discussed, the dark matter today is dominated by the axion in the region T γ,osc < T γ,c . The glueballs are dominant dark matter today in the region T γ,c < T γ,osc .
There are various astrophysical observations to constraint the mass of the glueball and axion dark matter.
We shortly summarize the bounds on it. When the glueball dominates dark matter, its self-interaction gives observable effects if the scattering rate is large enough to reach the isothermal density profile inside the halo. It can be also detectable from the merger of dark matter halos, because the glueballs will be slowed down during the collision, which leads to the offset between the dark matter and the collisionless components like stars. The crosssection of the glueball like self-interacting dark matter is bounded as ( [51] and references therein) In terms of the glueball mass, the glueball should have the mass greater than O(50)MeV if it is the dominant component of dark matter. The phenomenology of heavier glueball dark matter was studied in [52]. If the axion is the dominant dark matter component, there is also the lower bound on the axion mass due to its fuzziness. The de Broglie wavelength becomes astrophysical scale if m a is around 10 −22 eV, and suppresses the structure formation. The ultra-light axion dark matter can act like waves that are bound to or interact with each other by gravity inside the halo, which leads to the formation of solitonic cores and macroscopic quasiparticles moving around the center. These structures can have a great influence on the motion of stars. All these considerations give the strong constraint on the axion mass in the range m a 10 −22 − 10 −20 eV ( [53,54] and references therein). There is another constraint due to the observation of highly spinning black holes, even if the axion is not dark matter. That is because if the axion mass is close to the inverse of the size of the spinning black hole, a superradiance phenomenon occurs and parts of black hole's mass and spin are removed by the produced axion cloud. Current observations of the supermassive black hole with the mass of 10 6 − 10 7 M provide interesting constraints within the axion mass range 10 −20 − 10 −16 eV ( [55] and reference therein). However, it is sensitive to the axion self-interaction, the constraint is model dependent.
When the dark matter is mostly composed of the axions, Ω a h 2 0.11, and the freeze-out of the glueball happens after the BBN, the fraction of the dark glueball subcomponent dark matter becomes One can think that there is no strong constraint on the subcomponent glueball dark matter for f g 0.1. However, the evolution of the glueball dark matter after structures form may alter the cosmological history of the Universe around z = 7 − 15 as discussed in Sec. V.

IV. PERTURBATIONS
We now study the evolution of the cosmological perturbations for the axion and glueball dark matter. Both have non-trivial features compared to the CDM. For example, the transition of the nature of the axion from dark energy to dark matter modifies the early ISW effect of the cosmic microwave background [56]. The perturbation at scales smaller than the effective de Broglie wavelength is suppressed by its wave nature [3,57,58]. For glueballs, the number-changing self-interaction also disturbs the growth of the density perturbation of the scales which enter the horizon well before the freeze-out [11,59].
On one hand, in our set-up, the dark sector is decoupled from the visible sector and the origin of their abundance can be totally different from that of the SM particles. Let us provide a simple example for it. As the origin of Eq. (2.1), the axion can originate from the phase of a complex scalar field X. The corresponding matter Lagrangian at high scales is Here Q is the vector-like fermion charged under the dark gauge group. We have the anomalous global symmetry: After X gets a nonzero vacuum expectation value from the scalar potential, it can be decomposed as The axion φ is massless at perturbative level, and the radial scalar s gets a mass as m s = λ/2f a , which can be much smaller than f a for λ 1. The mass of dark fermion is M Q = yf a / √ 2. At scales µ < M Q , integrating out the heavy fermion yields the following Lagrangian The interaction between s (φ) and dark gluons is coming from the one-loop diagram mediated by the heavy fermion Q. Finally the radial scalar s is decoupled at scales µ m s , and we arrive at Eq. (2.1). Since the dark sector is decoupled from the visible sector, the initial amount of dark sector particles can be generated by a completely different mechanism than that of the Standard Model particles. For example, if the inflation Hubble parameter H I is given as m s H I M Q , the U (1) PQ symmetry is not restored, and s and a are all light particles during inflation. The stochastic fluctuations of the scalar field during inflation δs ∼ H I /2π can yield its initial abundance after inflation. From Eq. (4.4), these scalars will eventually decay mostly to axions and partially to dark gluons as Br(s → gg) ∼ (N g 2 /8π 2 ) 2 with a total decay rate Γ s ∼ m 3 s /8πf 2 a . Gluons will be quickly thermalized and form thermal bath with a temperature T g , while the axions are just redshifted. Although this is just one of the production mechanisms for dark sector, it gives a good motivation to study the isocurvature perturbation of dark matter from the initial dark gluon temperature fluctuation δT g,i . In this case, in addition to the adiabatic perturbation, there are two sources of the isocurvature perturbation. One is the fluctuation of the axion field during the inflation δφ i , and the other is the fluctuation of the gluon temperature δT g,i , which is induced by initial perturbation of the decaying scalar as Eq. Based on the evolution of the background dark axion and dark gluon/glueball, we solve the equations for the density perturbations focusing on the effect of isocurvature perturbation transfer and obtain the approximated solutions for the super-horizon modes (k = 0), in order to understand the parametric dependence more clearly.

A. Adiabatic perturbation
For the evolution of the multicomponent fields or fluids, the perturbations can be decomposed into the curvature (adiabatic) perturbation and the isocurvature (entropy) perturbations. The adiabatic perturbation is the modes perturbed along the direction of the background evolution, so that for any different species X and Y [60][61][62]. S XY is the relative entropy perturbation, whose name can be easily understood from thermodynamics. For an isolated species which satisfies the continuity equation ρ X + 3H(ρ X + p X ) = 0, the perturbation of the entropy density s X is given as δs X /s X = Hδρ X /ρ X , hence S XY = 3 δ ln(s X /s Y ). The adiabatic mode can be described by the evolution of the comoving curvature perturbation [63], The corresponding initial condition for the adiabatic mode is derived as where δ γ is for the photon fluid, the index i indicates the time at which the initial perturbation is defined. Here we use the axion potential and mass in Eq. (2.8).
From Eq. (2.18) and (2.20), we derive the solutions for the super-horizon modes in radiation-dominated era.
where m a (T g (τ )) ≡ dm a (T g (τ ))/dτ . As the scale factor becomes larger than a ci , the terms proportional to m a is rapidly vanishing. One can show that from the continuity equation for the coupled gluon-axion fluid, the detailed evolution of m a (T g ) around a = a ci does not lead to the different final result. Eq. (4.8) states that the adiabatic perturbation shares same form as δ X = (1 + w X )R i after the confining phase transition of the dark sector, and no history dependence happens, because S XY = 0 holds under the time evolution for the super-horizon modes. This is the characteristic feature of the adiabatic perturbation.

B. Isocurvature perturbation
The isocurvature perturbation is a mode perturbed along a direction orthogonal to the direction of background evolution. Taking S X as we can trace the evolution of the individual component of the isocurvature perturbation. At the linear perturbation level, the curvature perturbation cannot generate the isocurvature perturbations, and it is also conserved on the super-horizon scales [61]. Thus, for the evolution of isocurvature perturbations, we can safely take R i = 0, so that Φ i = 0, Ψ i = 0, and δρ tot,i = 0 as the initial condition at high temperatures, and solve the perturbation equations of the density contrast of X with the initial nonzero δ X,i The actual evolution of the density perturbation can be numerically calculated based on Eq. (2.25) and compared with the CMB and matter power spectrum. The form of the perturbation becomes particularly simple, if both the axion and glueballs becomes CDM-like well before the matter-radiation equality. In this case, the constraints on the isocurvature perturbation can be easily provided by comparing the analytic formula in super-horizon limit with the criteria of the Planck 2018 [1]. Therefore in this section, let us focus on this case.
The isocurvature perturbation of dark matter is expressed as where 'hat' denotes the Gaussian random variables satisfying δ a,iδg,i = 0, and with T aa T gg 1. The off-diagonal elements of the transfer matrix, T ag and T ga are determined by Eq. (2.25) as follows.

Induced by the initial displacement of the axion field
The evolution of the isocurvature perturbation induced by an initial density perturbation of the axion δ a,i = 2 δφ i /φ i can be described by the input value of δ a,i with the condition R i = 0 and the associated solutions from Eq. (2.25c). For T g T g,c , ∂V /∂T g ≈ 0 and In the axion dominated dark matter scenario (Ω a Ω g ), the dominant contribution is trivial: (δ DM ) iso δ a,i . In the opposite case (Ω a Ω g ), the relevant equation for the glueball isocurvature perturbation induced by that of the axion is where c 2 g ≡ (dp g /d ln T g )/(∂ρ g /∂ ln T g ). Note that the combination of 1 + 3Hc 2 g T g /T g is vanishing in the limit of ρ a /ρ g = 0 as This implies that the transfer matrix element T ga , which represents the effect of the initial axion isocurvature perturbation to that of the glueball dark matter, is of O(Ω 2 a /Ω 2 g ). It is the consequence of the conservation of the gluon/glueball entropy during the evolution. Therefore, the dominant contribution to the isocurvature perturbation of dark matter is just that from the axion subcomponent dark matter. In summary, for Ω g Ω a , (δ DM ) iso Ω a Ω DMδ a,i for Ω a Ω g . (4.16)

Induced by the initial fluctuation of the gluon temperature
For the initial fluctuation of the gluon temperature δ g,i ≈ 4 δT g,i /T g,i , the condition R i = 0 and δφ i = 0 for the perturbative variables give the following initial values In the glueball-dominated dark matter case, the contribution of the axion is suppressed by its energy density, so (δ DM ) iso δ g,i . In the opposite case (Ω g Ω a ), the effect of the gluon temperature fluctuation to the axion perturbation can be captured by Eq. (2.22). There are three stages of the axion evolution: (I) slow-rolling period (H m a (T g )), (II) confining phase transition with the saturation of the axion mass m a (T g ) = m a , (III) axion oscillating period (H m a ). For (I), the axion mass term is negligible and The solution becomes δφ φ η a 2(2η a + 4)(2η a + 5) After the confining phase transition, the perturbation of the axion in the periods (II), (III) obeys the equation of motion without δm a term, δφ + 2Hδφ + a 2 m 2 a δφ = 0. with the constant coefficients δφ ± . Matching Eq. (4.19) and Eq. (4.21) at a = a ci determines δφ ± and gives the solution for (II) and (III). In the period (III), δ a is given by where H ci is the Hubble rate at a = a ci . Note that the transfer matrix element T ag is suppressed by the factor of m 2 a /H 2 ci whenever the axion dominates the dark matter, but no further suppression happens. Therefore, for Ω a Ω g , (4.23) C. Bound on the isocurvature perturbation Since δφ i and δT g,i are independent random fluctuations, the power spectrum can be decomposed as P (k, z) = P RR (k, z) + X=a,g P II,X (k, z), (4.24) where X = a, g stand for the isocurvature perturbations induced by δ a,i and δ g,i , respectively. For the decomposition of the power spectrum as P (k, z) = (2π 2 /k 3 )P(k)T (k, z), we can match the primordial spectrum P(k) with the values we obtained for the superhorizon modes in the previous section. From the observations, the adiabatic mode is nearly scale-independent. For the isocurvature perturbations produced during the inflation, we can also naturally assume they are nearly scale-invariant, where k * is the pivot scale of the wave number, and parameters {A X , n X } are constants.
Observation of the CMB presents the upper bound on the isocurvature perturbation [1,[64][65][66]. The constraint is expressed by the bound on the isocurvature fraction β iso , which is defined by , (4.26) where P RR and P II are the power spectra defined in Eq. (4.25). For density perturbations, P RR = k 3 R 2 i /2π 2 and P II = X k 3 (δ DM ) 2 iso,X /2π 2 . We focus on the large scales to constrain the primordial perturbation from the CMB data. The constraint on β iso for the pivot scale is given by [1] β iso (k * = 0.002 Mpc −1 ) < 0.035. (4.27) We can compare this constraint with the values of β iso in our scenario. If the axion dominates dark matter, Ω g Ω a Ω DM (R(r, f a ) < 1 for Eq. (3.26)), the fraction β iso can be expressed as Because R(r, f a ) < 1, the term proportional to Ω g /Ω DM is always the dominant contribution in the second term of the RHS. In the opposite case, if the glueball is the main dark matter component, Ω a Ω g Ω DM (R(r, f a ) > 1), we have where . (4.31) Although the coupling between the axion and the gluon is key to the amount of axion dark matter, the contributions of the same coupling to the isocurvature perturbation is always subdominant. The perturbation is close to the sum of the independent elements, so it allows a large isocurvature perturbation of the subcomponent dark matter. The effect of such a large isocurvature perturbation is not clear yet. Since the glueballs are strongly self-interacting particles, it may provide non-trivial effects if the glueball is the subcomponent dark matter with a large isocurvature perturbation.
In the following section, we study the somewhat different aspect of the subcomponent glueball dark matter in the late time Universe.
In the standard mechanism on the formation and growth of black holes, SMBHs efficiently increase their mass by the accretion of baryonic material. However, the rate is limited, because the radiation pressure slows down the absorption. The maximal growth rate is captured by the Salpeter time based on the Eddington limit [67,68], where is the efficient factor, m p is the proton mass. If the seed black hole is generated at t i with a mass M seed , If the seed black hole is formed at z = 15, the maximal black hole mass becomes (2 − 6) × 10 4 M seed at z = 7. If the seed is formed at z = 30, its mass becomes (6 − 10) × 10 5 M seed . Therefore in order to explain the SMBHs with masses of O(10 9 M ) at z ∼ 7, a seed mass should be greater than (10 4 − 10 5 )M . This is quite challenging in the standard theory of black hole formation.
On the other hand, by solving the gravo-thermal fluid equations [14] and performing N -body simulation [15] with the assumption that the host halo is isolated, it is shown that such a heavy seed black hole could be generated from the gravo-thermal collapse of the subcomponent dark matter. Given the NFW density profile for the dominant DM component, the seed black hole is formed with the mass when the age of the Universe becomes t(z col ) = t(z i ) + ∆t col .
Here, M h is the mass of the host halo, z corresponds to the redshift. t(z i ) is the time when the virialized dark matter halo is isolated as we assume. ∆t col is the duration of the gravo-thermal collapse of the subcomponent dark matter for given initial conditions. β 1 and ∆t col are both calculated numerically. The fraction factor β 1 0.025/(ln(1 + c) − c/(1 + c)) in [14], where c is the concentration of the NFW profile (M h = 4πρ s r 3 s (ln(1 + c) − c/(1 + c))), and β 1 0.006 in [15]. By comparing the dark matter halo density profiles of two papers, we find that both results are well matched. The formation period ∆t col is estimated as the form where β 2 456 (480), p = 0 (2) in [14] ( [15]), and the apparent relaxation time of the subcomponent dark matter at t = t(z i ) is defined as = 0.28 Myr 10 cm 2 /g f g σ g /m g Note that ∆t col can be shorter than the age of the Universe for a given z, t(z) 550 Myr( 10 1+z ) 3/2 . Therefore, for the isolated halo with a mass M h = 10 12 M , f p+1 g σ g /m g (1 − 10) cm 2 /g, and f g 0.001 − 0.01 can explain the SMBH around z = 7. We illustrate the formation of the seed black hole and its growth history in Fig. 3 for the halo mass M h = 10 12 M .
In our scenario, the dark glueball dark matter provides a such strongly interacting subcomponent dark matter.  [14,15]. All information in red illustrates parameter space for a seed black hole (red dot). The seed black hole can be on the Eddington curve or on the shaded area in which the observations are explained by slower growth of the seed black hole. The time of collapse (z col ) and the mass of the seed black hole M seed are determined by model parameters {fg, σg/mg} or {mg, r}.
Since β 2 and p are directly estimated in N -body simulation, we take the result of [15] (β 2 = 480, p = 2) as the benchmark value. Then, the relevant combination of the model parameters is f 3 g σ g /m g , which is estimated as For the final expression, Eq. (3.35) is used. Because it is very sensitive to r, the ratio parameter is nearly predicted from the explanation of the SMBH at high z. The corresponding allowed range of the glueball mass is also provided as m g = O(0.05)MeV for f g = O(0.001). As to the parameters of the dominant component of dark matter, the axion, its decay constant is f a = O(10 16 GeV) and the axion mass becomes m a = O(10 −18 ) eV. This is safe from the current fuzzy dark matter constraints. Interestingly, this axion mass is also related with the supermassive black hole with the mass of M BH ∼ 10 7 M through superradiance as we discussed before. The axions can be efficiently generated from the spinning black hole by superradiant amplification. During the amplification, the axion also takes away the sizable amount of the black hole's angular momentum, which gives the contradiction to the observation [69]. However, if the selfinteraction among the axions is sizable, they will collapse before the axion cloud is saturated [70], and the loss of the angular momentum is limited. For m a ∼ 10 −18 eV, the GUT scale decay constant provides a sizable axion self-interaction to trigger bosenovae. Therefore, the constraint may not be applied directly.
Several simplications are used in the previous discussion. Let us discuss possible caveats and alternative history of the seed black hole formation. The host halo mass is taken as 10 12 M . This is because the halo mass is expected to be greater than O(10 3 ) times the mass of its SMBH [71,72]. In N -body simulations [73][74][75], the comoving number density of the cold dark matter halos with M h ≥ 10 12 M is evaluated as (10 −5 −10 −6 )(Mpc) −3 at z = 7. Thus, the halo is also heavy enough to coincide with the fact that observations of SMBH around z = 7 are rare.
However, since we consider the formation of the seed black hole at higher redshifts (z > 7), the existence of such (isolated) heavy halo is questionable. If we extrapolate the halo mass function obtained by the N -body simulation [75], the comoving number density of the halos with M h ≥ 10 12 M becomes (10 −8 − 10 −9 )(Mpc) −3 at z = 10, and 10 −15 (Mpc) −3 at z = 15. In this context, the issue of formation of heavy seed black holes is just transferred to the problem of supermassive halo formation at high redshifts.
On one hand, based on N -body simulations, we can define M h (z) at a given z in such a way that the comoving number density of the halos with their masses greater than M h (z) is given by 10 −6 (Mpc) −3 . Then, M h (z) is evaluated as 10 12 M at z = 7, 10 11 M at z = 10, and 10 10 M at z = 15. It is more natural to think the possibility that when the seed black hole is formed, the mass of the host halo is smaller than 10 12 M , although it is still one of the heaviest halos at z i . These heaviest halos get bigger and bigger by mergers with nearby smaller halos or by accretion of the gases. The actual merger history is quite complex, but the heaviest halo is likely to remain the heaviest. In this sense, we consider M h (z) as the evolution of the host halo mass, and estimate the growth rate Γ h (z) as . (5.10) The last equality holds numerically for 7 z 15. The black hole growth rate by the accretion of baryons is much greater than the halo growth rate. However, the halo mass is still hierarchically larger than the black hole mass during the evolution. The another important feature is that in terms of the halo mass, the relaxation time defined by Eq. (1 + z) The concentration parameter c also depends on the halo mass and the redshift. The recent N -body simulation [76] calculates the concentration parameter c(M h , z) as the function of M h and z in a wide range of M h Therefore, we expect that the gravo-thermal collapse for the SMBH would not be triggered if m g keV. If m g is much larger than O(keV), the three-to-two interaction is not effective before the gravo-thermal collapse happens. The gravo-thermal collapse begins to accelerate after ∆t col . During the collapse, the diffusion of the dark matter mass is inefficient, and the glueballs concentrate their mass of O(M seed ) around the center by increasing the core density and its temperature [77,78]. Then, the number changing interaction becomes gradually important. It is not clear how it affects the last stage of the gravo-thermal collapse (the formation of the seed black hole). This is because the temperature increase rate caused by gravo-thermal collapse is not known yet for such a high mass density of the core. We leave it for future work.
There is also the lower bound on the glueball mass from the cosmological evolution. If the glueball is light enough, it becomes a warm or hot dark matter, so that its speed around z = 7 − 15 is greater than the escape velocity of the halo. This means that the subcomponent dark matter is not clustered, and cannot provide a good initial condition. After the dark glueball freeze-out, its velocity scales as 1/a. The corresponding redshifted glueball velocity at a given z is v g (z) 10 −3 1 + z 16 r 0.001 3 2 100 eV m g 5 4 , (5.18) if the freeze-out happens before the epoch of matterradiation equality, and v g (z) 10 −3 1 + z 16 r 0.001 4 3 100 eV m g 10 9 , (5.19) if the freeze-out happens in the dark matter dominated era. In order to explain the SMBH formation, this value should be hierarchically smaller than the virial velocity v s ∼ 10 −3 during the period z = 7 − 15. In our scenario, r is nearly fixed as 0.005. See Eq. (5.9). This implies the lower bound on m g as 100 eV.

VI. DISCUSSION
We have studied the cosmological evolution of dark light scalars, whose masses and interactions originate from the approximate global symmetry and the nonperturbative dynamics of the hidden gauge symmetry. One is the feebly interacting dark axion, and the other is the strongly interacting dark glueball. Both can be dark matter if they are light enough. The equations of motion are derived and evaluated to identify the dark matter abundance and the perturbation evolution induced by the coupling between the axion and the dark gluon. We also explore the possibility that the subcomponent glueball dark matter contributes to the formation of the supermassive black hole at redshift z ∼ 7.
Although we have dealt with the problems as closely as possible, there are still many questions that have not been covered by this paper. What would be the observable consequences of the first-order confining phase transition? In our discussion, we ignore gravitational wave productions during the confining phase transition, because it is weakly first-order unless N is very large. However, if the phase transition happens around the recombination era, it may leave a footprint on the CMB. What is the exact form of the axion scalar potential and the effect of self-interactions? The scalar potential of the axion is not a simple cosine form, and a multi-branch structure may provide the nontrivial effects if the axion is produced around the spinning supermassive black hole by superradiant amplification. What is the correct period of the gravo-thermal collapse when the fraction of the subcomponent dark matter is small enough? So far, there is no intensive study on the gravo-thermal collapse of the subcomponent dark matter for such a small fraction. The empirical form of the collapse time scale ∆t col should be confirmed for f g 10% and higher scattering cross-sections. What is the effect of the number changing interactions of the glueball dark matter for the final stage of the black hole formation? During the gravo-thermal collapse, one may think of the possibility that the defining phase transition occurs, because of the large density of the glueball dark matter inside the core. It would be very interesting to study the implication of such a microscopic nature of the dark matter for the final formation of the black hole.