# $O(10)M_\odot$ primordial black holes and string axion dark matter

Keisuke Inomata, Masahiro Kawasaki, Kyohei Mukaida, Yuichiro Tada, Tsutomu T. Yanagida

► To cite this version:

Keisuke Inomata, Masahiro Kawasaki, Kyohei Mukaida, Yuichiro Tada, Tsutomu T. Yanagida. *$O(10)M_\odot$ primordial black holes and string axion dark matter*. Physical Review D, 2017, 96 (12), pp.123527. 10.1103/PhysRevD.96.123527. hal-01703688

HAL Id: hal-01703688

https://hal.science/hal-01703688v1

Submitted on 5 Oct 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
O(10)M☉ primordial black holes and string axion dark matter

Keisuke Inomata,<fnref n="1" />,<fnref n="2" /> Masahiro Kawasaki,<fnref n="1" />,<fnref n="2" /> Kyohei Mukaida,<fnref n="2" /> Yuichiro Tada,<fnref n="1" />,<fnref n="2" />,³,⁴ and Tsutomu T. Yanagida<fnref n="2" />

¹ICRR, University of Tokyo, Kashiwa 277-8582, Japan

²Kavli IPMU (WPI), UTIAS, University of Tokyo, Kashiwa 277-8583, Japan

³Institut d'Astrophysique de Paris, UMR-7095 du CNRS, Université Pierre et Marie Curie,
98 bis bd Arago, 75014 Paris, France

⁴Sorbonne Universités, Institut Lagrange de Paris, 98 bis bd Arago, 75014 Paris, France
(Received 10 October 2017; published 26 December 2017)

LIGO-Virgo collaboration has found black holes as heavy as $M \sim 30M_☉$ through the detections of the gravitational waves emitted during their mergers. Primordial black holes (PBHs) produced by inflation could be an origin of such BHs. While it is tempting to presume that these PBHs constitute all dark matter (DM), there exist a number of constraints for PBHs with $O(10)M_☉$ which contradict with the idea of PBHs as all DM. Also, it is known that weakly interacting massive particle (WIMP) that is a common DM candidate is almost impossible to coexist with PBHs. These observations motivate us to pursue another candidate of DM. In this paper, we assume that the string axion solving the strong CP problem makes up all DM, and discuss the coexistence of string axion DM and inflationary PBHs for LIGO events.

DOI: 10.1103/PhysRevD.96.123527

# I. INTRODUCTION

Axion appearing in string theory [1–3] is a plausible solution to the strong CP problem, which may provide a platform to discuss why there exists an approximate global symmetry with an extreme accuracy in the low-energy theory. The decay constant of axion is around the grand unified theory scale, $f_a \sim 10^{16}$ GeV, in many string axion models. Although a mild tuning of the misalignment angle is required not to overclose the Universe, the coherently oscillating axion can account for the present dark matter (DM) abundance.<fnref n="1" /> Since axion is light, they may develop super-horizon fluctuations during inflation. To avoid severe bounds on the isocurvature perturbation [6], the Hubble parameter during inflation has to be low enough $H_{inf} \lesssim 10^9$ GeV($f_a/10^{16}$ GeV)$^{0.405}$.

However, low-scale inflation has a serious drawback: namely an extreme fine-tuning of the initial condition is required (See, e.g., [7]). Interestingly, from the viewpoint of the string landscape, such an initial condition can be selected dynamically. In the string landscape, there are many different vacua with different vacuum energies, and inflations may occur in each vacuum [8–11]. A preceding inflation takes place with a higher vacuum energy, ends by tunneling or violating the slow-roll condition, and the Universe falls into a vacuum with a smaller energy. After successive such processes, eventually, the low-scale inflation responsible for our observable patches of the Universe takes place. To be concrete, let us consider hilltop inflation that is a prominent example of low-scale inflation. If the inflaton for hilltop inflation acquires a positive

Hubble-induced mass term, the inflaton can stay at the top of the flat potential during those pre-inflations. As a result, the infamous initial condition problem is solved dynamically [12], which completes a consistent cosmological scenario for the string axion DM.

Recently, the LIGO-Virgo collaboration has announced the detections of three gravitational wave (GW) events, GW150914 [13], GW151226 [14], and GW170104 [15]. In addition to these three events, there is one GW candidate called LVT151012, whose signal-to-noise ratio is somewhat smaller than those of the other three GWs [14]. These GWs are generated by mergers of two black holes (BHs). Their masses are $36.2_{-3.8}^{+5.2} M_☉$ and $29.1_{-4.4}^{+3.7} M_☉$ in GW150914, $14.2_{-3.7}^{+8.3} M_☉$ and $7.5_{-2.3}^{+2.3} M_☉$ in GW151226, $31.2_{-6.0}^{+8.4} M_☉$ and $19.4_{-5.9}^{+5.3} M_☉$ in GW170104, and $23_{-6}^{+18} M_☉$ and $13_{-5}^{+4} M_☉$ in LVT151012. One can see that the three BHs out of the eight BHs are as heavy as $30 M_☉$. However, in the usual metallicity ($Z \sim Z_☉$), stellar origin BHs may not be $30 M_☉$ due to mass loss by stellar wind [16–18]. Currently, researchers are actively searching for the origin of these BHs.

Primordial black hole (PBH) is one of the candidates of such BHs [19–23].<fnref n="2" /> It can be formed by the collapse of overdense Hubble patch in the early radiation-dominated universe [25–27], which is completely different from the standard formation mechanism of stellar BHs with super novae. That is why PBH can have various masses such as $30 M_☉$. Also, it is clear that inflation can be an origin of such large density perturbations [28–30]. To realize PBHs

<footnote n="1">¹See Refs. [4,5] for discussion on the entropy production to avoid such a tuning.</footnote>

<footnote n="2">²In a low-metallicity environment, stellar origin BHs with $30M_☉$ can be produced [16–18]. This kind of BH is another candidate of the $30M_☉$ BHs [24].</footnote>
however, the curvature perturbations on small scales should be strongly amplified roughly by $10^4$ times compared to those on the cosmic microwave background (CMB) scale, though such a rapid change of the perturbation amplitude is generally difficult in the slow-roll single-field inflation [31]. This is because the tilt of the power spectrum of the curvature perturbations are determined by the slow-roll parameters, which are small during slow-roll inflation. On the other hand, if we have another inflation during the last 50–60 e-folds as supported by the viewpoint of the string landscape, the second inflation can be free from the COBE normalization. Therefore, although fine-tuning of the parameters is required as we show in Sec. IV, such an extreme modulation of perturbations can be achieved easily compared to ordinary single-field inflations [29,32–36].

Moreover, the axion DM and PBHs can coexist, while the weakly interacting massive particle (WIMP) cannot [37,38]. Though the minimal setup might be the case where PBHs for the LIGO events account for the present DM simultaneously, there are many constraints on this PBH mass range which strongly disfavor PBHs as all DM. Hence, we focus on the possibility that PBHs and other DM candidate coexist. Since PBHs themselves and/or primordial perturbations large enough for abundant PBHs cause compact DM clumps [39–43], annihilations of DM in the present Universe can be dramatically enhanced, which makes the mixture of PBHs and WIMP DM difficult. In this sense, the axion is suitable because its interaction is feeble.

In this paper, we assume that $\mathcal{O}(10)M_\odot$ BHs detected by LIGO-Virgo collaboration are PBHs, and construct a concrete scenario where the string axion DM and PBHs for the LIGO events coexist. The organization of this paper is as follows. In Sec. II, we review the string axion DM and the isocurvature perturbations produced by the fluctuations of an axion field. We also discuss the properties of the inflation model compatible with the string axion DM. In Sec. III, we discuss PBHs for the LIGO events. The DM candidate in the presence of $\mathcal{O}(10)M_\odot$ PBHs is also discussed. In Sec. IV, we construct the scenario consistent with the string axion DM. We take the double inflation model as a concrete example and show the concrete parameter consistent with the observations. Section V is devoted to the conclusion.

## II. STRING AXION DM AND ISOCURVATURE PERTURBATION

In this section, we summarize basic cosmological properties of the string axion. Throughout this paper, we adopt the following normalization of the axion interaction:

$$ \mathcal{L}_{aF\tilde{F}} = \frac{g^2}{16\pi^2 f_a} a \mathrm{Tr} F_{\mu\nu} \tilde{F}^{\mu\nu}. \quad (1) $$

For many string axion models, the axion decay constant tends to be large, and the value is around the GUT scale [1–3]:

$$ f_a \sim 10^{16} \text{ GeV}. \quad (2) $$

We take it as a fiducial value.

### A. String axion as DM

Suppose that the moduli field (saxion) is stabilized throughout the evolution of our Universe. Around the QCD phase transition, the string axion acquires a mass term via the QCD instanton, and its potential becomes $\sim m_a^2(T)a^2$ with $m_a(T)$ being the temperature dependent axion mass. After this mass term exceeds the Hubble parameter, it starts to oscillate and behaves as cold dark matter (CDM). In this case, the current energy density is determined by the initial amplitude, namely the summation of the initial misalignment and the superhorizon fluctuations during inflation [44]:

$$ \Omega_a h^2 \approx 10^4 \left( \frac{f_a}{10^{16} \text{ GeV}} \right)^{1.19} (\theta^2 + \sigma_\theta^2), \quad (3) $$

where $\theta$ is the misalignment angle and $\sigma_\theta^2$ is the variance of the fluctuations given by$^3$

$$ \sigma_\theta^2 \approx \left( \frac{H_{\mathrm{inf}}}{2\pi f_a} \right)^2. \quad (4) $$

The present DM density is $\Omega_{\mathrm{DM}}h^2 \approx 0.12$ [45]. By requiring $\Omega_a \le \Omega_{\mathrm{DM}}$ and using Eq. (3), we can set an upper bound on the misalignment angle [5,44],

$$ \theta \lesssim 3.4 \times 10^{-3} \left( \frac{f_a}{10^{16} \text{ GeV}} \right)^{-0.595}, \quad (5) $$

and the Hubble parameter during inflation,

$$ H_{\mathrm{inf}} \lesssim 2 \times 10^{14} \text{ GeV} \left( \frac{f_a}{10^{16} \text{ GeV}} \right)^{0.405}. \quad (6) $$

As we show below, the constraint on the isocurvature perturbation puts a more stringent bound on the Hubble parameter during inflation.

### B. Isocurvature perturbation

The axion condensate starts to oscillate when its mass exceeds the Hubble parameter. Thus, for a homogeneous initial field value of axion, its density perturbation is necessarily adiabatic. However, since the mass of axion

<footnote>$^3$Precisely speaking, there is an additional logarithmic factor, $\ln(k_{\mathrm{max}}L)$ with $k_{\mathrm{max}}$ and $L$ being a comoving wave number entering the horizon at the end of inflation and the infrared cut off respectively. Though, parametrically, its value is $\mathcal{O}(10)$, the precise value is not so important for our purpose because we expect $\theta \gg \sigma_\theta$ if the string axion dominates the present DM density, as we will see in the subsequent discussion.</footnote>
is light, the axion acquires fluctuations during inflation and the fluctuations yield differences in the initial condition for each Hubble patch at the onset of its oscillation. As a result, the axion also develops the isocurvature perturbation. The power spectrum of the axion isocurvature perturbation can be estimated as [44]

$$P_{\text{iso}} \approx 4 \left( \frac{\Omega_a}{\Omega_{\text{DM}}} \right)^2 \frac{(H_{\text{inf}}/2\pi f_a)^2}{\theta^2 + \sigma_{\theta}^2}. \quad (7)$$

By analyzing the CMB spectrum, the Planck collaboration sets an upper bound on the isocurvature perturbation:

$$\beta_{\text{iso}} \equiv \frac{P_{\text{iso}}}{P_{\text{ad}} + P_{\text{iso}}} < 0.038 \quad (8)$$

at 95% CL for the scale of $k = 0.05 \text{ Mpc}^{-1}$ [46]. Here $P_{\text{ad}}(\approx 2.1 \times 10^{-9})$ is the power spectrum of the adiabatic perturbation. For the case where the string axion constitutes all the DM, we can derive an upper bound on the Hubble parameter during inflation from Eqs. (7) and (8) as [5]

$$H_{\text{inf}} \lesssim 10^9 \text{ GeV} \left( \frac{f_a}{10^{16} \text{ GeV}} \right)^{0.405} . \quad (9)$$

Note that, in the case where the string axion with $f_a \sim 10^{16}$ GeV makes up all the DM, we expect $\theta \gg \sigma_\theta$, and $P_{\text{iso}} \approx 4(H_{\text{inf}}/2\pi f_a\theta)^2$. This also means that the isocurvature constraint, Eq. (9), sets a severer bound than that of Eq. (6).

### C. Hilltop inflation after preinflations

We have seen that the inflation scale has to be low ($H_{\text{inf}} \lesssim 10^9$ GeV) to avoid the severe constraints on the isocurvature perturbation for the string axion. Hilltop (new) inflation [47–49] is one of the prominent examples of the low-scale inflation and thus we consider it in the following.

The major drawback of hilltop inflation is its initial condition problem. If we assume that any inflationary phase does not precede hilltop inflation, then there seems to be no reason to expect that the inflaton is smoothened over the horizon and its expectation value lies near a local minimum of its potential. However, once we realize that there could be multiple inflations before the observed one, this problem can be solved dynamically. Also, the picture of multiple pre-inflations may be natural in the context of the string landscape [8–11].

To make our discussion concrete, let us consider a preinflation right before the observed one. An inflaton $\varphi$ which is responsible for the hilltop inflation may acquire a positive Hubble-induced mass term during the preinflation via $V_{\text{pr}}\varphi^2/M_{\text{Pl}}^2$ with $V_{\text{pr}}$ being a potential for the preinflation. If this is the case, $\varphi$ is stabilized at the local minimum throughout the preinflation, and hence the initial condition is set dynamically. A concrete realization can be found in Ref. [12]. There, it was assumed that a single

hilltop inflation after the preinflation explains the whole e-folds. In general, there might be two (or more) low-scale inflations that account for the required e-folds in total [29,32–36]. This consideration opens up a possibility to generate PBHs as we will see in the following sections.

### III. PRIMORDIAL BLACK HOLES AND THE LIGO EVENTS

In this section, we discuss a possibility to explain the LIGO GW events by mergers of PBHs and the required property for the primordial density perturbation to have such PBHs. We also discuss DM candidates which can coexist with PBHs.

#### A. LIGO GW events

The LIGO-Virgo collaboration has detected three GW events and one candidate from mergers of BH-binary. Though the mass distribution is still unclear, those BHs have masses of $O(10)M_{\odot}$. They have also estimated the merger rate as $12–213 \text{ Gpc}^{-3} \text{yr}^{-1}$ [15]. The question is an amount of PBHs required to explain these events. This is addressed in Ref. [21]; the fraction of PBHs has to be around $O(10^{-3})–O(10^{-2})$ to reproduce the merger rate.<fnref n="4" /> Unfortunately, the estimated event rate might suggest that those PBHs do not constitute all the DM. However, at the same time, the required amount of PBHs is still consistent with several observational constraints [52–61].

#### B. PBH formation

Here we briefly review the formation of PBHs by the primordial density perturbation and discuss the suitable property of the primordial density perturbation to account for the LIGO events.

Basically, PBHs are formed by the gravitational collapse of an over-dense region. According to simple analysis by Carr [27], if a density contrast is larger than $\delta > 1/3$ at its horizon reentry, the over-dense region overcomes the radiation pressure and PBHs are formed. We use this threshold value $\delta_c = 1/3$ as a fiducial value in the following. The PBH mass is roughly estimated as the horizon mass at the horizon reentry of the perturbations which are larger than the threshold. We can estimate the relation between the scale of the perturbations and PBH mass as

$$ \begin{aligned} M(k) &= \gamma \rho \frac{4\pi H^{-3}}{3} \Big|_{k=aH} \approx \frac{\gamma M_{\text{eq}}}{\sqrt{2}} \left(\frac{g_{*,\text{eq}}}{g_*}\right)^{\frac{1}{6}} \left(\frac{k_{\text{eq}}}{k}\right)^2 \\ &\approx M_{\odot} \left(\frac{\gamma}{0.2}\right) \left(\frac{g_*}{10.75}\right)^{-\frac{1}{6}} \left(\frac{k}{1.9 \times 10^6 \text{ Mpc}^{-1}}\right)^{-2} \end{aligned} \quad (10) $$

<footnote n="4">⁴There are some papers that discuss the constraints on the $O(10)M_{\odot}$ PBH abundance from the merger rate and the constraints are $\Omega_{\text{PBH}}/\Omega_{\text{DM}} < O(10^{-2})$ [50,51].</footnote>
$$ \approx M_{\odot} \left(\frac{\gamma}{0.2}\right) \left(\frac{g_{*}}{10.75}\right)^{-\frac{1}{6}} \left(\frac{f}{2.9 \times 10^{-9} \text{ Hz}}\right)^{-2}, \quad (11) $$

where $\gamma$ is a numerical factor which depends on the properties of gravitational collapse and $M_{eq}$ is the horizon mass at matter-radiation equality time and $g_*$ and $g_{*,eq}$ are the degrees of freedom at PBH formation and matter-radiation equality time, where $g_* \approx 10.75$ at $\mathcal{O}(10)M_\odot$ PBH production. We adopt the simple analysis value $\gamma = (1/3)^{3/2} \approx 0.2$ as a fiducial value [27].

Assuming that the density perturbations follow the Gaussian distribution, the production rate of PBHs with $M_{\text{PBH}} = M$ can be written as

$$ \beta(M) = \int_{\delta_c} \frac{d\delta}{\sqrt{2\pi}\sigma^2(M)} e^{-\frac{\delta^2}{2\sigma^2(M)}} \approx \frac{1}{\sqrt{2\pi}\delta_c/\sigma(M)} e^{-\frac{\delta^2}{2\sigma^2(M)}}. \quad (12) $$

$\sigma^2(M)$ is the variance of density contrast with the smooth- ing scale $M(k)$ which is defined as [62]

$$ \sigma^2(M(k)) = \int d\ln q W^2(qk^{-1}) \frac{16}{81} (qk^{-1})^4 P_{\zeta}(q), \quad (13) $$

where $P_{\zeta}(k)$ is the power spectrum of the primordial curvature perturbation. $W(x)$ is a window function and throughout this paper we use the Gaussian window function $W(x) = e^{-x^2/2}$.

Once PBHs are formed, the PBH energy density behaves as nonrelativistic matter and decays as $\rho_{\text{PBH}} \propto a^{-3}$. On the other hand, the background energy density behaves as radiation $\rho \propto a^{-4}$ until the matter-radiation equality time. Therefore the ratio of the PBH energy density to the background energy density grows as $\rho_{\text{PBH}}/\rho \propto a$. Using the PBH production rate $\beta(M)$, we can derive PBH-DM ratio over logarithmic mass interval $d\ln M$ as

$$ \begin{aligned} & \frac{\Omega_{\text{PBH}}(M)}{\Omega_{\text{DM}}} \\ & \approx \left. \frac{\rho_{\text{PBH}}(M)}{\rho_m} \right|_{\text{eq}} \frac{\Omega_m}{\Omega_{\text{DM}}} = \left( \frac{T_M}{T_{\text{eq}}} \frac{\Omega_m}{\Omega_{\text{DM}}} \right) \gamma \beta(M) \quad (14) \\ & \approx \left( \frac{\beta(M)}{1.84 \times 10^{-8}} \right) \left( \frac{\gamma}{0.2} \right)^{\frac{3}{2}} \left( \frac{10.75}{g_*(T_M)} \right)^{\frac{1}{4}} \left( \frac{0.12}{\Omega_{\text{DM}}h^2} \right) \left( \frac{M}{M_\odot} \right)^{-\frac{1}{2}}, \quad (15) \end{aligned} $$

where we define $\Omega_{\text{PBH}}(M) \equiv d\Omega_{\text{PBH}}/d\ln M$ and $\rho_{\text{PBH}}(M) \equiv d\rho_{\text{PBH}}/d\ln M$. $\Omega_m$ and $\Omega_{\text{DM}}$ are the energy density parameter of total matter (baryon + CDM) and CDM. $T_M$ and $T_{\text{eq}}$ are the temperatures of the Universe at the horizon reentry of the perturbations and the matter-radiation equality time. We adopt $\Omega_{\text{DM}}h^2 = 0.12$ [45]. The total PBH-DM ratio is given by

$$ \Omega_{\text{PBH,tot}} = \int d\ln M \Omega_{\text{PBH}}(M). \quad (16) $$

To sum up, we need a large curvature perturbation, $P_{\zeta}(k) \sim 10^{-2}$, at a scale of $k \sim 10^6$ Mpc$^{-1}$ to reproduce the LIGO events by mergers of PBHs. Since the curvature perturbation at the large scale, $k \lesssim 1$ Mpc$^{-1}$, should be as small as $\sim 10^{-9}$ not to conflict with the CMB observation, we need to break the (almost) scale-invariance of the spectrum. In addition, there are other constraints on such a large curvature perturbation at the scale, $k \sim 10^6$ Mpc$^{-1}$, as pointed out in Ref. [63–65]; (i) CMB spectrum distortion and (ii) induced GWs by the second order effect. Basically, the peak of the power spectrum of the curvature perturbations should be sufficiently sharp$^5$ to avoid them and one must explicitly check that the model is well consistent with those constraints to investigate a successful scenario.

### C. DM candidates in the presence of PBHs

Throughout this paper, we focus on the string axion in a low-scale inflation scenario. Thus, our primary DM candidate is axion. Nevertheless, it is instructive to compare the axion DM with other DM candidates in the presence of PBHs.

The fascinating possibility would be the case where PBHs for the LIGO events can constitute the whole DM simultaneously. Unfortunately, it seems to be difficult to achieve full DM because of many observational constraints at this mass range [52–61]. Hence, we need to consider other DM candidates.

The weakly interacting massive particle (WIMP) is a major particle DM candidate. However, it is known that WIMP is difficult to coexist with PBHs. This is because a large density perturbation promotes the formation of DM halos with a steep profile, which is called ultracompact minihalos (UCMHs) [39–43]. As a result, the DM density at the center of the UCMH becomes so dense that DM annihilations at the present Universe are dramatically enhanced. Thus, we can constrain such scenarios by the observation of gamma-rays from DM annihilations. The upper bound on the annihilation cross section is known to be several orders of magnitude smaller than its typical value for WIMP, $\langle\sigma v\rangle \sim 10^{-26}$ cm$^3$ s$^{-1}$ [38]. Moreover, even if PBHs are not produced by large primordial perturbations, the existence of PBHs itself leads to the formation of UCMHs with PBHs as their cores. Therefore the PBH abundance is constrained as $\Omega_{\text{PBH}} \lesssim 10^{-4}$ (for $m_{\text{DM}} \sim 100$ GeV) for a vast range in PBH mass with the thermal relic DM [37]. These facts indicate the difficulty of the coexistence of PBHs and WIMP DMs.

<footnote>$^5$One can avoid such constraints also by an enhanced non-Gaussianity [66,67].</footnote>
One important constraint on the coexistence of the string axion and BHs comes from the superradiance effect [68–70]. We briefly mention it here. If a Compton wavelength of axion is close to the size of a rotating BH, then such a BH immediately loses its spin by filling bounded Bohr orbits with a huge number of axions due to the bosonic nature of axions [71–73]. As a result, observations of rotating BHs can constrain axion models [74,75]. Parametrically, we expect that a given rotating BH mass can constrain $m_a \sim 1/r_g = 4\pi M_{Pl}^2/M_{BH} \sim 10^{-11}$ eV($M_☉/M_{BH}$), or equivalently $f_a \sim 10^{18}$ GeV($M_{BH}/M_☉$). One can see that ordinary stellar BHs with $M_{BH} \sim M_☉$ cannot probe the axion decay constant much below the Planck scale. Hence, $f_a \sim 10^{16}$ GeV in our scenario is still viable.

## IV. DOUBLE INFLATION AFTER PREINFLATIONS

As mentioned at the end of Sec. II, if we have two (or more) inflations which explain the required e-folds, the primordial curvature perturbation can be large enough for a sizable amount of PBHs by breaking the scale-invariance and generate PBHs. In this paper, we consider the double inflation model [29] as a concrete example.<fnref n="6" /> In this section, first, we search for the desired parameter regions semi-analytically, taking account of two conditions. One is that the resonance enhancement of perturbations are small enough for a linear analysis. The other one is that the first inflation should correctly reproduce the result of the Planck observations on the CMB scale. After showing the desired parameter regions, we will show the fully numerical result with one desired parameter set in the last subsection.

### A. Inflation scenario

The double inflation model has two stages of inflation. During the intermediate phase between the first inflation and the second inflation, this model can generate the large perturbations which produce PBHs. Note that the large-scale perturbations observed by the Planck satellite are generated by the first inflation and the small-scale perturbations are generated by the second inflation.

As discussed in Sec. II, the string axion DM requires the low energy scale of the first inflation as

$$H_{\text{inf,first}} \lesssim 10^9 \text{ GeV}, \quad (17)$$

because the isocurvature perturbations generated by the first inflation are severely constrained by CMB observations. In the following, we assume that the first and second inflations are low-scale hilltop inflations to satisfy Eq. (17).

We take the following inflaton potentials:

$$V(\varphi, \chi) = V_1(\varphi) + V_2(\chi) + V_{\text{stb}}(\varphi, \chi) + C, \quad (18)$$

$$V_1(\varphi) = v_1^4 \left( -\frac{\varphi}{\varphi_l} - \frac{\varphi^2}{\varphi_q^2} + \left( 1 - \frac{\varphi^{n_1}}{\varphi_{\min}^{n_1}} \right)^2 \right), \quad (19)$$

$$V_2(\chi) = v_2^4 \left( -\frac{\chi}{\chi_l} - \frac{\chi^2}{\chi_q^2} + \left( 1 - \frac{\chi^{n_2}}{\chi_{\min}^{n_2}} \right)^2 \right), \quad (20)$$

$$V_{\text{stb}}(\varphi, \chi) = \frac{1}{2} c_{\text{pot}} V_1(\varphi) \chi^2, \quad (21)$$

where $M_{Pl}$ is set to be unity and $\varphi$ is the inflaton of the first inflation and $\chi$ is the inflaton of the second inflation. $C$ is a compensate term which makes $V(\varphi_{min}, \chi_{min}) = 0$.<fnref n="7" /> We take $n_1, n_2 \ge 3$.<fnref n="8" /> Generally speaking, we also expect Planck-suppressed order correction to kinetic terms,

$$\mathcal{L}_{\text{kin}} = -\frac{1}{2} \left( 1 - \frac{c_{\text{kin}}}{2} \chi^2 \right) \partial_{\mu} \varphi \partial^{\mu} \varphi - \frac{1}{2} \partial_{\mu} \chi \partial^{\mu} \chi + \dots \quad (22)$$

The expected inflation dynamics is described as follows. First, the first inflation takes place while $\varphi$ slowly rolls down the potential and the energy scale of the first inflation is $\sim v_1^4$. During the first inflation, the other inflaton $\chi$ is stabilized at $\chi \approx v_2^4/(c_{\text{pot}}v_1^4\chi_l)$ because of the stabilization term, $\frac{1}{2}c_{\text{pot}}V_1(\varphi)\chi^2$. After a while, the slow-roll of $\varphi$ ends and $\varphi$ starts to oscillate around the potential minimum, $\varphi_{\min}$. During this time, $\varphi$ possibly passes many times through the tachyonic region, where the mass of $\varphi$ becomes tachyonic, and nonadiabatic region, where the adiabatic condition is violated, and the perturbations of $\varphi$ can grow so large that the perturbations become nonlinear. We discuss this issue again in the next subsection and Appendix A. After the amplitude of the oscillation becomes small due to the Hubble friction, the energy density deposited in $\varphi$ becomes smaller than the second inflation energy scale $\sim v_2^4$ and then the second inflation starts. By then, the stabilization of $\chi$ becomes so weak that $\chi$ can roll down the potential. After the slow-roll of $\chi$ ends, $\chi$ oscillates around the minimum and decays to other particles at least via Planck-suppressed operators. Around minimum, the mass of $\chi$ is given by

<footnote n="7">⁷Strictly speaking, the true potential minimum is not $V(\varphi_{min}, \chi_{min})$. However, the deviation of the true potential minimum from $V(\varphi_{min}, \chi_{min})$ is so small and we have checked that the deviation does not affect our results.</footnote>

<footnote n="8">⁸If we take $n_1, n_2 = 2$, we must take the value of inflaton at the potential minimum as $\varphi_{min} > M_{Pl}$ or $\chi_{min} > M_{Pl}$ to realize slow-roll inflation. If an inflaton becomes larger than Planck scale, the inflaton potential can be modified by the other Planck suppressed operators we do not write explicitly and it can be difficult to follow the dynamics of the inflaton. Therefore we restrict our discussion to the case of $n_1, n_2 \ge 3$.</footnote>

<footnote n="6">⁶While our model causes double inflation with two fields, some inflation models cause double inflation with single field [76,77]. In [76,77], they discuss PBHs with their mass around $O(10^{20})$ g.</footnote>
$$m_{\chi} = \sqrt{2}n_2 \frac{v_2^2}{\chi_{\min}}. \qquad (23)$$

Then the corresponding reheating temperature can be evaluated as

$$T_R \approx 0.1m_{\chi}^{3/2} \approx 0.2n_2^{3/2} \frac{v_2^3}{\chi_{\min}^{3/2}}, \qquad (24)$$

if the reheating is achieved by some Planck-suppressed operators. We use this value in the following analysis.

Let us estimate the e-folds which corresponds to the peak scale of curvature perturbations. We define the pivot scale as $k_{\text{pivot}} = 0.05 \text{ Mpc}^{-1}$ and the peak scale as $k_{\text{peak}} = 5 \times 10^5 \text{ Mpc}^{-1}$ (see Fig. 2). Then the difference between the corresponding e-folds is calculated as

$$N_{\text{pivot}} - N_{\text{peak}} = \log\left(\frac{k_{\text{peak}}}{k_{\text{pivot}}}\right) \approx 16. \qquad (25)$$

Note that there is a difference between the scale at the first inflation end and the peak scale as discussed in Appendix B. The difference between $N_{\text{peak}}$ and $N_{\text{first inf,end}}$ (= e-folds at the end of the first inflation) is given by

$$\frac{H_{\text{inf,first}}}{H_{\text{inf,second}}} \approx e^{3(N_{\text{peak}} - N_{\text{first inf,end}})}, \qquad (26)$$

where $H_{\text{inf,first}}$ and $H_{\text{inf,second}}$ are the energy scales of the first inflation and the second inflation. Here we define $N_{\text{first,pivot}}$ as $N_{\text{first,pivot}} = N_{\text{pivot}} - N_{\text{first inf,end}}$. If we take

$H_{\text{inf,first}}/H_{\text{inf,second}} = 400$, we derive $N_{\text{peak}} - N_{\text{first inf,end}} \approx 2$ and $N_{\text{first,pivot}} = 18$. In the following, we use $H_{\text{inf,first}}/H_{\text{inf,second}} = 400$ and $N_{\text{first,pivot}} = 18$ as fiducial values.

## B. Resonance

Some inflation models have the preheating phases with some resonances, which follow the inflation phases. In general, hilltop inflation has a preheating phase driven by the tachyonic mass and nonadiabaticity after its inflation. Throughout this paper, we focus on the preheating phase which occurs after the first inflation because we must follow the evolution of the perturbation which exits the horizon during the second inflation. Roughly speaking, the resonance becomes stronger for smaller $\varphi_{\min}$ and larger $H_{\text{inf,first}}$ because the ratio of the oscillation timescale to the Hubble friction timescale ($H_{\text{inf,first}}/m_{\varphi} = \sqrt{3/2}\varphi_{\min}/n$) becomes smaller for smaller $\varphi_{\min}$ and the seed perturbation $H_{\text{inf,first}}/2\pi$ becomes larger for higher $H_{\text{inf,first}}$. In Fig. 1, we show the parameter regions where the perturbations become nonlinear, which is estimated by numerical calculations. We discuss the behavior of the resonances further in Appendix A. If the resonances are strong and the perturbations become nonlinear, it becomes complicated to calculate the evolution of perturbations because the evolution depends on the higher order perturbations and backreactions of perturbations to the homogeneous backgrounds become non-negligible. In the double inflation model we consider here, there are the interaction terms between inflatons $\varphi$ and $\chi$, such as $\frac{1}{2}c_{\text{pot}}V_1(\varphi)\chi^2$ and $\frac{1}{4}c_{\text{kin}}\chi^2\partial_{\mu}\varphi\partial^{\mu}\varphi$. Hence, we must follow the dynamics of

<figure><img src="image_2.png" /><figcaption>FIG. 1. The blue shaded regions show the parameter regions where both $n_s$ and $dn_s/d\ln k$ are consistent with observations. In each point of this figure, we calculate the value of $H_{\text{inf,first}}$ with which $A_s$ is consistent with the observation. The lower left white region in the left figure shows the parameter region of $N_{\text{first,pivot}} < 18$. The light black shaded regions show the parameter regions in which perturbations become nonlinear due to the resonances (see also Appendix A for details). We take $\varphi_q = \sqrt{200}$ and $N_{\text{first,pivot}} = 18$ as a fiducial value and we take $n = 3$ in the left figure and $n = 4$ in the right figure.</figcaption></figure>
the inflaton $\varphi$ during its oscillation phase to calculate the perturbations. In this paper, to avoid the complexities, we search for the parameters in which the perturbations remain linear until horizon exit. The study of the nonlinear regime is left for future works.

### C. CMB constraints

Now, we return to the double inflation model and discuss the constraints on the large-scale perturbations. From the Planck results, the cosmological parameters are determined in 95% CL as [45]:

$$10^9 A_s = 2.142 \pm 0.098, \qquad (27)$$

$$n_s = 0.9667 \pm 0.0080, \qquad (28)$$

$$\frac{dn_s}{d \ln k} = -0.0085 \pm 0.0152, \qquad (29)$$

$$r_{0.002} < 0.113, \qquad (30)$$

where $A_s$ is the amplitude of the adiabatic scalar power spectrum, $n_s$ is the scalar spectral index, $dn_s/d \ln k$ is the running of scalar spectral index, and $r_{0.002}$ is the tensor-to-scalar ratio.

Since the first inflation generates the curvature perturbations on the pivot scale, we must check whether the first inflation can be consistent with Eqs. (27)–(30). Thus we focus on the first inflation whose potential is given by Eq. (19).<fnref n="9" />

The slow-roll parameters, $\varepsilon$, $\eta$ and $\zeta$, are defined as

$$\varepsilon \equiv \frac{1}{2} \left( \frac{V_1'}{V_1} \right)^2, \qquad (31)$$

$$\eta \equiv \frac{V_1''}{V_1'}, \qquad (32)$$

$$\zeta \equiv \frac{V_1' V_1'''}{V_1^2}, \qquad (33)$$

where a prime denotes a derivative with respect to $\varphi$. The cosmological parameters can be written with the slow-roll parameters as

$$A_s \equiv \frac{H_{\text{inf,first}}^2}{8\pi^2\varepsilon}, \qquad (34)$$

$$n_s \approx 1 - 6\varepsilon - 2\eta, \qquad (35)$$

$$\frac{dn_s}{d \ln k} \approx -24\varepsilon^2 + 16\varepsilon\eta - 2\zeta, \qquad (36)$$

$$r_{0.002} \approx 16\varepsilon. \qquad (37)$$

We can calculate the $e$-folds from the end of the first inflation with the following equation:

$$N(\varphi) = \int_{\varphi_e}^{\varphi} \frac{V}{V'} d\varphi, \qquad (38)$$

where $\varphi_e = (\varphi_{\text{min}}^n/(2n_1(n_1-1)))^{\frac{1}{n_1-2}}$. According to the discussion in the previous subsection about resonance, we define $\varphi_{\text{first,pivot}}$ as the value satisfying $N(\varphi_{\text{first,pivot}}) = 18$. If we determine $\varphi_l$, $\varphi_q$ and $\varphi_{\text{min}}$, we can calculate $\varphi_{\text{first,pivot}}$ numerically. Then we can substitute it into Eqs. (31)–(37) and derive $A_s$, $n_s$, $dn_s/d \ln k$ and $r_{0.002}$. Since the relation $\varepsilon \ll \eta$ is generally satisfied in hilltop inflation models, Eq. (30) is always satisfied in our concrete inflation model and we do not consider Eq. (30) explicitly in the following. Here, let us stress that $n_s$ and $dn_s/d \ln k$ are independent of $H_{\text{inf,first}}$. This is because the slow-roll parameters, $\varepsilon$, $\eta$ and $\zeta$, are independent of $H_{\text{inf,first}}$ as

$$\varepsilon \approx \frac{1}{2} \left( -\frac{1}{\varphi_l} - 2 \frac{\varphi}{\varphi_q^2} - \frac{2n}{\varphi_{\text{min}}} \left( \frac{\varphi}{\varphi_{\text{min}}} \right)^{n_1-1} \right)^2, \qquad (39)$$

$$\eta \approx -\frac{2}{\varphi_q^2} - \frac{2n_1(n_1-1)}{\varphi_{\text{min}}^2} \left( \frac{\varphi}{\varphi_{\text{min}}} \right)^{n_1-2}, \qquad (40)$$

$$\zeta \approx \left( -\frac{1}{\varphi_l} - 2 \frac{\varphi}{\varphi_q^2} - \frac{2n}{\varphi_{\text{min}}} \left( \frac{\varphi}{\varphi_{\text{min}}} \right)^{n_1-1} \right) \times \left( -2 \frac{2n(n-1)(n-2)}{\varphi_{\text{min}}^3} \left( \frac{\varphi}{\varphi_{\text{min}}} \right)^{n_1-3} \right). \qquad (41)$$

To satisfy the slow-roll condition, $\eta$ must be sufficiently smaller than one. For simplicity, we take $\varphi_q = \sqrt{200}$ in the following. Note that while $n_s$ and $dn_s/d \ln k$ depend on only $\varphi_l$, $\varphi_q$ and $\varphi_{\text{min}}$, $A_s$ depends also on $H_{\text{first,inf}}$ in addition to those parameters.

In Fig. 1, the blue shaded regions show the parameter regions of $\varphi_l$ and $\varphi_{\text{min}}$ where the corresponding $n_s$ and $dn_s/d \ln k$ are consistent with the observation. In each point of this figure, we calculate the value of $H_{\text{inf,first}}$ with which $A_s$ is consistent with the observation. The parameter regions which cause the strong resonance are shown by the light black shaded regions. One can see that in order for perturbations to avoid the strong resonance and for $n_s$ and $dn_s/d \ln k$ to be consistent with the observational values, $H_{\text{inf,first}}$ must be

<footnote n="9">⁹The relation between cosmological parameters and the parameters of this inflation potential has also been studied in [78,79].</footnote>
$$H_{\text{inf,first}} \gtrsim \begin{cases} 10^8 \text{ GeV} & \text{for } n=3, \\ 3.2 \times 10^{10} \text{ GeV} & \text{for } n=4. \end{cases} \quad (42)$$

If we take $n > 4$, the lower bound of $H_{\text{inf,first}}$ becomes higher because the perturbations in the case of the larger $n$ are accompanied by the stronger resonance.

## D. Concrete example

In this subsection, let us show the full numerical result for the following parameters:

$$\begin{align}
n_1 &= 3, & v_1 &= 2.667 \times 10^{-5} \quad (H_{\text{inf,first}} = 10^9 \text{ GeV}), \nonumber \\
\varphi_l &= 7.088 \times 10^5, & \varphi_q &= \sqrt{200}, & \varphi_{\min} &= 0.1169, \nonumber \\
n_2 &= 3, & v_2 &= 1.334 \times 10^{-6} \quad (H_{\text{inf,second}} = 2.5 \times 10^6 \text{ GeV}), \nonumber \\
\chi_l &= 1.196 \times 10^{11}, & \chi_q &= \sqrt{10}, & \chi_{\min} &= 5.223 \times 10^{-3}, \nonumber \\
c_{\text{pot}} &= 0.675, & c_{\text{kin}} &= -0.671, & C &= 3.465 \times 10^{-23}, \tag{43}
\end{align}$$

as a fiducial example in the desired parameter regions obtained in the previous subsection.

In Fig. 2, the resultant linear perturbations are indicated with the constraints by the $\mu$-distortion and BBN. Though we chose the parameters not to yield the resonance, the power spectrum shows a sharp peak on $k \sim 5 \times 10^5 \text{ Mpc}^{-1}$. This is because we apply the amplification mechanism proposed in Ref. [80], that is, $c_{\text{pot}}$ and $c_{\text{kin}}$ are chosen so that the effective Hubble-induced mass of $\chi$ during the $\varphi$-oscillation phase, $m_{\chi}^2 \approx \frac{3}{2}(c_{\text{pot}} + c_{\text{kin}})H^2$, becomes sufficiently small. Then the perturbation of $\chi$ can avoid the damping, which leads to the effective amplification of

<figure><img src="image_3.png" /><figcaption>FIG. 2. The scalar power spectrum for parameters given in Eq. (43). Orange shaded regions are excluded by the current constraint on $\mu$-distortion, $|\mu| < 9 \times 10^{-5}$ [81] and the effect on n-p ratio during big-bang nucleosynthesis [82] (There are other constraints on the power spectrum of curvature perturbations around $10^4-10^5 \text{ Mpc}^{-1}$ [83,84] and their constraints are similar to [82].) The black dotted line represents a future constraint by $\mu$-distortion with the PIXIE [85], $|\mu| < 10^{-8}$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3. The PBH mass spectrum for parameters given in Eq. (43). Red shaded regions are excluded by extragalactic gamma-rays from Hawking radiation (EGY) [86], femtolensing of known gamma-ray bursts (Femto) [87], white dwarfs existing in our local galaxy (WD) [88], microlensing search with Subaru Hyper Suprime-Cam (HSC) [89], Kepler microlensing (Kepler) [90], EROS/MACHO microlensing (EROS/MACHO) [91], dynamical heating of ultra faint dwarf galaxies (UFD) [53], and accretion constraints from CMB (CMB) [52] (Recently, Poulin et al. have revisited the accretion constraints from CMB taking account of nonspherical accretion and got $\Omega_{\text{PBH}}/\Omega_{\text{DM}} \lesssim O(10^{-2})$ around $O(10)M_{\odot}$ [92].) See also [23] for a recent summary of observational constraints on PBHs and [58–61] for the constraints on the extended PBH mass spectrum.</figcaption></figure>

$\delta\chi$. See Appendix B in Ref. [80] for a more detailed description. Such a sharp peak is required to avoid the pulsar timing array (PTA) constraints as we will see below.

Using equations described in Sec. III, we calculate the expected PBH abundance from the obtained curvature perturbations and show its result in Fig. 3.<fnref n="10" /> From this, one can see that there is a peak at $M_{\text{PBH}} = 30 M_{\odot}$ and the height is around $\Omega_{\text{PBH}}/\Omega_{\text{DM}} \approx O(10^{-3})$. According to [21], the expected merger rate in this PBH spectrum is consistent with the merger rate estimated by LIGO-Virgo collaboration, $12-213 \text{ Gpc}^{-3} \text{ yr}^{-1}$ [15].

As pointed out in Refs. [64,65], the density (scalar) perturbations large enough to realize abundant PBHs can generate GWs (tensor perturbations) due to the second order effect, and the frequency of such secondary GWs for $O(10)M_{\odot}$ PBHs lies around the sensitivity region of the PTA experiments ~1 nHz. Therefore one has to carefully check the consistency with these experiments. Following the instruction of Ref. [80], we calculate the current abundance of the secondary GWs and show the result in Fig. 4 with the current and future observational constraints.

<footnote n="10">¹⁰In Fig. 3, we plot the monochromatic observational constraints by the red shaded regions. Strictly speaking, since our PBH mass spectrum is an extended mass function, the constraints can be tighter as discussed in [58–61]. To take into account the constraints on the extended mass spectrum, we have used the analysis described in [58,61] and confirmed that our PBH mass spectrum is consistent with the observational constraints.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 4. The induced GW spectrum for parameters given in Eq. (43). Green shaded regions are excluded by PTA observations with EPTA [93], PPTA [94], and NANOGrav [95]. A black dotted line shows the prospect of the SKA sensitivity [96,97].</figcaption></figure>

Thanks to the sharpness of our power spectrum of the
curvature perturbations, the resultant GWs successfully
evade the current PTA constraints.

V. CONCLUSIONS

In this paper, we have discussed $O(10)M_\odot$ PBHs in the
presence of string axion DM.

The axion naturally appears in the context of string
theory and can be DM through its coherent oscillation.
Since the string axion is light, its perturbations which exit
the horizon during inflation behave as isocurvature pertur-
bations constrained by the CMB observations. The
presence of string axion DM puts the severe constraint
on the Hubble parameter during the inflation as $H_{inf} \lesssim$
$10^9$ GeV($f_a/10^{16}$ GeV)$^{0.405}$ and therefore the energy scale
of inflation must be low.

We have constructed the scenario where the string axion DM and PBHs for the LIGO events coexist, taking the double (hilltop + hilltop) inflation model. We have shown that this inflation model can satisfy the following conditions at the same time.

(i) The energy scale of the inflation model must be low enough to be consistent with the constraint on the isocurvature perturbations.

(ii) The large-scale perturbations observed by CMB must be consistent with the observational values of $A_s$, $n_s$, $d n_s / d \ln k$ and $r_{0.002}$.

(iii) The perturbations which make $O(10)M_\odot$ PBHs must be large enough for the PBH-DM ratio to be $\Omega_{PBH}/\Omega_{DM} \approx O(10^{-3}) - O(10^{-2})$ and sharp enough to be consistent with the constraints from PTA, $\mu$-distortion, and BBN observations.

(iv) The perturbations must remain linear until the second inflation ends to be calculated with the linear analysis.

Note that the last condition does not necessarily mean that
the nonlinear perturbations cannot satisfy the other

conditions. The case with the nonlinear perturbations is
left for future works.

Finally, let us discuss several outcomes of our scenario.

First, our model predicts the sharp peak around $k \sim 10^6$ Mpc$^{-1}$. This peak is constrained from the large scale ($\mu$-distortion, BBN) and small scale (PTA). In other words, our model can be tested by future observations such as PIXIE ($\mu$-distortion) and SKA (PTA).

Second, let us mention the axion miniclusters [98]. It has
been suggested that axion can form a dense clump by,
e.g., large isocurvature perturbations associated with the
breakdown of $U(1)_{PQ}$ symmetry after inflation [98–100].
Although we have considered the scenario where $U(1)_{PQ}$
symmetry is broken before inflation, our scenario expects
the large adiabatic perturbations around $k \sim 10^6$ Mpc$^{-1}$
and hence the axion miniclusters might be produced.<fnref n="11" />
Since the perturbations are smaller than that in the scenario
where $U(1)_{PQ}$ is broken after inflation (but still large), we
can expect that the produced axion miniclusters in our
scenario are sparser [102]. As discussed in [103], the sparse
axion miniclusters tend to get disrupted by the encounters
with the stars in the Galactic disk and the fragments of
the destructed miniclusters form tidal streams. Although
there are still many uncertainties about this phenomenon, it
is possible that the tidal streams can be detected in the
future.

ACKNOWLEDGMENTS

This work is supported by Grant-in-Aid for Scientific
Research from the Ministry of Education, Science, Sports,
and Culture (MEXT), Japan, No. 15H05889 (M. K.),
No. 17K05434 (M. K.), No. 17H01131 (M. K.),
No. 26104009 (T. T. Y.), No. 26287039 (T. T. Y.) and
No. 16H02176 (T. T. Y.), World Premier International
Research Center Initiative (WPI Initiative), MEXT, Japan
(K. I., M. K., K. M., Y. T., and T. T. Y.), JSPS Research
Fellowships for Young Scientists (Y. T., and K. M.), Région
Île-de-France (Y. T.), and Advanced Leading Graduate
Course for Photon Science (K. I.). T. T. Y. is a
Hamamatsu Professor.

APPENDIX A: RESONANCE

In this Appendix, we discuss the resonances which occur
after hilltop inflation.

For simplicity, we consider the toy inflation model with
the following potential:

$$
V(\phi) = v^4 \left( 1 - \frac{\phi^n}{\phi_{\min}^n} \right)^2 . \qquad (A1)
$$

<footnote n="11">¹¹Tkachev has remarked the possibility that the dense axion clumps are produced by the large primordial perturbations in Ref. [101].</footnote>
<figure><img src="image_6.png" /><figcaption>FIG. 5. The schematic image of the inflation potential given by Eq. (A1). This figure is a modified version of Fig. 1 in [104].</figcaption></figure>

If we take $v$ and $\phi_{\min}$ as $v_1$ and $\phi_{\min}$, this potential is a good approximation of Eq. (19) during the oscillation of $\phi$ which follows the first inflation. In this inflation model, after inflation, $\phi$ starts to oscillate around $\phi_{\min}$. We define $\phi_e$ as the value at the end of inflation and $\phi_m$ as the value which satisfies $V''(\phi_m) = 0$, where the prime denotes a derivative with respect to $\phi$, e.g., $V''(\phi) = \partial^2 V(\phi)/\partial\phi^2$. Figure 5 shows the schematic image of the hilltop inflation potential. We can derive the value of $\phi_e$ and $\phi_m$ as

$$ \phi_e = \left( \frac{\phi_{\min}^n}{2n(n-1)} \right)^{\frac{1}{n-2}}, \qquad (A2) $$

$$ \phi_m = \left( \frac{n-1}{2n-1} \right)^{\frac{1}{n}} \phi_{\min}, \qquad (A3) $$

where we have used the equality $|V''(\phi_e)/V| = 1$ to derive $\phi_e$.

This toy inflation model goes through the following steps. When $\phi < \phi_e$, $\phi$ slowly rolls down its potential and inflation occurs with $H_{\inf} = v^2/\sqrt{3}$. After $\phi$ passes through $\phi_e$, the slow-roll inflation ends and $\phi$ starts to oscillate around $\phi_{\min}$. At first, the amplitude of oscillations is large enough so that $\phi$ can return to the tachyonic region, $\phi < \phi_m$. During this period, the perturbations of $\phi$ grow rapidly due to the tachyonic resonance. Later, the amplitude becomes small due to the expansion of the Universe and $\phi$ cannot return to the tachyonic region. Although the tachyonic resonance stops, perturbations still grow due to the nonadiabatic resonance. Adiabatic condition, $\omega^2 \gg \dot{\omega}$ ($\omega^2 \equiv V''(\phi)$), fails violently when $\phi \sim \phi_m$ because $V''(\phi) \sim 0$ at $\phi \sim \phi_m$. Finally, the amplitude of $\phi$ becomes so small that $\phi$ does not pass through the nonadiabatic region and the perturbations decay due to Hubble friction.

Figure 6 shows the time evolution of power spectrum of $\delta\phi$ divided by $\phi_{\min}^2$ at the resonance peak scale, $P_{\delta\phi}(k = k_{\text{reso,peak}})/\phi_{\min}^2$, with $n = 4$, $H_{\inf} = 10^9$ GeV,

<figure><img src="image_7.png" /><figcaption>FIG. 6. The time evolution of $P_{\delta\phi}(k = k_{\text{peak}})/\phi_{\min}^2$ with $n = 4$, $H_{\inf} = 10^9$ GeV, and $\phi_{\min} = 0.109$. We use the scale factor as time variables, where $a_{\text{inf,end}}$ is the scale factor at the inflation end.</figcaption></figure>

and $\phi_{\min} = 0.109$. $P_{\delta\phi}(k = k_{\text{reso,peak}})/\phi_{\min}^2$ roughly corresponds to the squared ratio of the perturbation to the background field value. If $P_{\delta\phi}(k = k_{\text{reso,peak}})/\phi_{\min}^2 > 1$, we regard the perturbation as nonlinear. From Fig. 6, we can see that the perturbation grows until $\log(a/a_{\text{inf,end}}) \sim 2$ due to the tachyonic and nonadiabatic resonance and after the growth ($\log(a/a_{\text{inf,end}}) \ge 2$), the perturbation decays due to Hubble friction. Here, we define the conformal time when a perturbation becomes a maximum as $\eta_{\text{peak}}$, which corresponds to the conformal time at $\log(a/a_{\text{inf,end}}) \sim 2$ in Fig. 6.

Tachyonic resonance has been already discussed analytically in [104]. Following Ref. [104], we can estimate the peak scale of the tachyonic resonance as

$$ \frac{k_{\text{reso,peak}}}{a_{\text{inf,end}}} \approx (\phi_{\min})^{\frac{n-2}{2n}} m_{\phi} = \sqrt{2} n v^2 (\phi_{\min})^{-\frac{n+2}{2n}}, \qquad (A4) $$

where $k_{\text{reso,peak}}$ is the comoving wave number at the peak, $m_{\phi}$ is the mass of $\phi$ around $\phi_{\min}$ defined as $m_{\phi} = \sqrt{2}nv^2/\phi_{\min}$ and $a_{\text{inf,end}}$ is the scale factor at the inflation end. Since the resonances occur so rapidly and the peak scale is determined at the beginning of the resonances, $a_{\text{inf,end}}$ appears in Eq. (A4). We confirm that our numerical results are consistent with Eq. (A4). Note that the resonance peak scale is different from the peak scale of the curvature perturbation which produces PBHs in our model ($k_{\text{reso,peak}} \neq k_{\text{peak}}$).

Figure 7 shows the $\phi_{\min}$ dependence of $P_{\delta\phi}(k = k_{\text{reso,peak}}, \eta = \eta_{\text{peak}})/\phi_{\min}^2$. From this figure, we can see that the decrease of $\phi_{\min}$ or the increase of $n$ leads to the increase of the peak height. This results can be understood as follows. The peak height depends on the relation between the timescale of the oscillation and that of Hubble friction. The timescale of the oscillation corresponds to $m_{\phi}^{-1} = \phi_{\min}/(\sqrt{2}nv^2)$. On the other hand, the
<figure><img src="image_8.png" /><figcaption>FIG. 7. The $\phi_{\min}$ dependence of $P_{\delta\phi}(k = k_{\text{reso,peak}}, \eta = \eta_{\text{peak}})/\phi_{\min}^2$. We take $H_{\text{inf}} = 10^9$ GeV. The blue solid (dotted) line shows $P_{\delta\phi}(k = k_{\text{reso,peak}}, \eta = \eta_{\text{peak}})/\phi_{\min}^2$ in the case of $n=3$ ($n=4$). The gray shaded region shows $P_{\delta\phi}(k = k_{\text{reso,peak}}, \eta = \eta_{\text{peak}})/\phi_{\min}^2 > 1$, where perturbations become non-linear.</figcaption></figure>

timescale of Hubble friction is related to the Hubble parameter as $H_{\text{inf}}^{-1} \approx 1/(\sqrt{3}v^2)$. The important factor is the ratio of $m_{\phi}^{-1}$ to $H_{\text{inf}}^{-1}$. If we take a smaller value of $\phi_{\min}$ or a larger value of $n$, only $m_{\phi}^{-1}$ becomes smaller (but $H_{\text{inf}}^{-1}$ does not change). Then the inflaton can pass through the resonance region more times before its amplitude becomes small due to Hubble friction and the peak height becomes larger.

Finally, let us mention the $H_{\text{inf}}$ dependence of the peak height. To understand it analytically, we neglect the metric perturbations<fnref n="12" /> and then the equation of motion for $\delta\phi$ which corresponds to the peak scale can be written as

$$ 
\begin{aligned}
\ddot{\delta\phi} + 3H\dot{\delta\phi} + \frac{k_{\text{reso,peak}}^2}{a_{\text{inf,end}}^2}\delta\phi + V''(\phi)\delta\phi &= 0 \\
\Rightarrow \dot{\delta\phi} + 3H\dot{\delta\phi} + 2n^2v^4(\phi_{\min})^{-\frac{n+2}{n}}\delta\phi \\
&\quad + v^4\frac{\partial^2}{\partial\phi^2}\left(1-\left(\frac{\phi}{\phi_{\min}}\right)^n\right)^2\delta\phi = 0 \\
\Rightarrow \frac{d^2}{d\tau^2}\delta\phi + 3\tilde{H}\frac{d}{d\tau}\delta\phi + 2n^2(\phi_{\min})^{-\frac{n+2}{n}}\delta\phi \\
&\quad + \frac{\partial^2}{\partial\phi^2}\left(1-\left(\frac{\phi}{\phi_{\min}}\right)^n\right)^2\delta\phi = 0,
\end{aligned}
\quad (A5) $$

where $d/d\tau \equiv d/d(v^2t)$ and $\tilde{H} \equiv \frac{1}{a}\frac{da}{d\tau}$. From Eq. (A5), we can see that $H_{\text{inf}}(= v^2/\sqrt{3})$ determines only the time

scale of the perturbation evolution. Note that even though the growth rate of perturbations due to the resonances are independent of $H_{\text{inf}}$, the value of $\delta\phi$ before the resonances depends on $H_{\text{inf}}$ as $\delta\phi \approx H_{\text{inf}}/2\pi$. Therefore, we can estimate the $H_{\text{inf}}$ dependence as $P_{\delta\phi}(k = k_{\text{reso,peak}}, \eta = \eta_{\text{peak}})/\phi_{\min}^2 \propto H_{\text{inf}}^2$. This estimation is consistent with our numerical calculation which includes the metric perturbations. Using the result of Fig. 7 and this assumption about $H_{\text{inf}}$-dependence, we illustrate the model parameters with which the perturbations become nonlinear by the light black shaded regions in Fig. 1. We avoid these regions for a concrete example in this paper.

## APPENDIX B: RELATION BETWEEN THE PEAK SCALE AND THE FIRST INFLATION END SCALE

In this Appendix, we derive Eq. (26). We define $a_{\text{peak}}$ as the scale factor which satisfies $a_{\text{peak}}H_{\text{inf,first}} = k_{\text{peak}}$, $a_{\text{first,e}}$ as the scale factor at the first inflation end and $a_{\text{second,b}}$ as the scale factor at the second inflation beginning. The scale factor dependence of $aH$ is given by

$$ aH \propto \begin{cases} a & \text{during the first inflation,} \\ a^{-1/2} & \text{during the oscillation of inflaton } \phi. \end{cases} \quad (B1) $$

Taking account of the fact that $H$ is proportional to $a^{-3/2}$ during the oscillation phase, we get the following equation:

$$ \log\left(\frac{H_{\text{inf,first}}}{H_{\text{inf,second}}}\right) = \frac{3}{2}\log\left(\frac{a_{\text{second,b}}}{a_{\text{first,e}}}\right). \quad (B2) $$

Following the discussion of Appendix B in Ref. [80], we can estimate as $k_{\text{peak}} = a_{\text{second,b}}H_{\text{inf,second}}(= a_{\text{peak}}H_{\text{inf,first}})$. Using this estimation and Eq. (B1), we get the following equation:

$$ 
\begin{aligned}
\log\left(\frac{a_{\text{second,b}}}{a_{\text{first,e}}}\right) &= 2\log\left(\frac{a_{\text{first,e}}H_{\text{inf,first}}}{a_{\text{second,b}}H_{\text{inf,second}}}\right) \\
&= 2\log\left(\frac{a_{\text{first,e}}}{a_{\text{peak}}}\right).
\end{aligned}
\quad (B3) $$

From Eqs. (B2) and (B3), we can derive Eq. (26) as

$$ \log\left(\frac{H_{\text{inf,first}}}{H_{\text{inf,second}}}\right) = 3\log\left(\frac{a_{\text{first,e}}}{a_{\text{peak}}}\right) \quad (B4) $$

$$ \Rightarrow \frac{H_{\text{inf,first}}}{H_{\text{inf,second}}} = e^{3(N_{\text{peak}} - N_{\text{first end}})}. \quad (B5) $$

<footnote n="12">¹²We have numerically checked that the metric perturbations have little effect on the dynamics of $\delta\phi$ during the resonances.</footnote>
[1] J. P. Conlon, J. High Energy Phys. 05 (2006) 078.

[2] P. Svrcek and E. Witten, J. High Energy Phys. 06 (2006) 051.

[3] K. Choi and K. S. Jeong, J. High Energy Phys. 01 (2007) 103.

[4] M. Kawasaki, T. Moroi, and T. Yanagida, Phys. Lett. B **383**, 313 (1996).

[5] M. Kawasaki, T. T. Yanagida, and N. Yokozaki, Phys. Lett. B **753**, 389 (2016).

[6] M. Kawasaki and T. Yanagida, Prog. Theor. Phys. **97**, 809 (1997).

[7] R. Brandenberger, Int. J. Mod. Phys. D **26**, 1740002 (2017).

[8] S. Kachru, R. Kallosh, A. D. Linde, and S. P. Trivedi, Phys. Rev. D **68**, 046005 (2003).

[9] B. Freivogel and L. Susskind, Phys. Rev. D **70**, 126007 (2004).

[10] L. Susskind, arXiv:hep-th/0302219.

[11] R. Bousso and J. Polchinski, J. High Energy Phys. 06 (2000) 006.

[12] K. I. Izawa, M. Kawasaki, and T. Yanagida, Phys. Lett. B **411**, 249 (1997).

[13] B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. **116**, 061102 (2016).

[14] B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. **116**, 241103 (2016).

[15] B. P. Abbott et al. (VIRGO, LIGO Scientific), Phys. Rev. Lett. **118**, 221101 (2017).

[16] B. P. Abbott et al. (Virgo, LIGO Scientific), Astrophys. J. **818**, L22 (2016).

[17] K. Belczynski, T. Bulik, C. L. Fryer, A. Ruiter, J. S. Vink, and J. R. Hurley, Astrophys. J. **714**, 1217 (2010).

[18] M. Spera, M. Mapelli, and A. Bressan, Mon. Not. R. Astron. Soc. **451**, 4086 (2015).

[19] S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, Phys. Rev. Lett. **116**, 201301 (2016).

[20] S. Clesse and J. García-Bellido, Phys. Dark Univ. **15**, 142 (2017).

[21] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Phys. Rev. Lett. **117**, 061101 (2016).

[22] Yu. N. Eroshenko, arXiv:1604.04932.

[23] B. Carr, F. Kuhnel, and M. Sandstad, Phys. Rev. D **94**, 083504 (2016).

[24] T. Kinugawa, K. Inayoshi, K. Hotokezaka, D. Nakauchi, and T. Nakamura, Mon. Not. R. Astron. Soc. **442**, 2963 (2014).

[25] S. Hawking, Mon. Not. R. Astron. Soc. **152**, 75 (1971).

[26] B. J. Carr and S. W. Hawking, Mon. Not. R. Astron. Soc. **168**, 399 (1974).

[27] B. J. Carr, Astrophys. J. **201**, 1 (1975).

[28] J. Garcia-Bellido, A. D. Linde, and D. Wands, Phys. Rev. D **54**, 6040 (1996).

[29] M. Kawasaki, N. Sugiyama, and T. Yanagida, Phys. Rev. D **57**, 6050 (1998).

[30] J. Yokoyama, Phys. Rev. D **58**, 083510 (1998).

[31] H. Motohashi and W. Hu, Phys. Rev. D **96**, 063503 (2017).

[32] M. Kawasaki and T. Yanagida, Phys. Rev. D **59**, 043512 (1999).

[33] P. H. Frampton, M. Kawasaki, F. Takahashi, and T. T. Yanagida, J. Cosmol. Astropart. Phys. 04 (2010) 023.

[34] M. Kawasaki, A. Kusenko, and T. T. Yanagida, Phys. Lett. B **711**, 1 (2012).

[35] M. Kawasaki, K. Mukaida, and T. T. Yanagida, Phys. Rev. D **94**, 063509 (2016).

[36] M. Kawasaki, A. Kusenko, Y. Tada, and T. T. Yanagida, Phys. Rev. D **94**, 083523 (2016).

[37] B. C. Lacki and J. F. Beacom, Astrophys. J. **720**, L67 (2010).

[38] T. Bringmann, P. Scott, and Y. Akrami, Phys. Rev. D **85**, 125027 (2012).

[39] J. A. Fillmore and P. Goldreich, Astrophys. J. **281**, 1 (1984).

[40] E. Bertschinger, Astrophys. J. Suppl. Ser. **58**, 39 (1985).

[41] M. Ricotti and A. Gould, Astrophys. J. **707**, 979 (2009).

[42] A. D. Ludlow, J. F. Navarro, V. Springel, M. Vogelsberger, J. Wang, S. D. White, A. Jenkins, and C. S. Frenk, Mon. Not. R. Astron. Soc. **406**, 137 (2010).

[43] M. Vogelsberger, S. D. White, R. Mohayaee, and V. Springel, Mon. Not. R. Astron. Soc. **400**, 2174 (2009).

[44] M. Kawasaki and K. Nakayama, Annu. Rev. Nucl. Part. Sci. **63**, 69 (2013).

[45] P. A. R. Ade et al. (Planck), Astron. Astrophys. **594**, A13 (2016).

[46] P. A. R. Ade et al. (Planck), Astron. Astrophys. **594**, A20 (2016).

[47] L. Boubekeur and D. Lyth, J. Cosmol. Astropart. Phys. 07 (2005) 010.

[48] A. D. Linde, Phys. Lett. B **108B**, 389 (1982).

[49] A. Albrecht and P. J. Steinhardt, Phys. Rev. Lett. **48**, 1220 (1982).

[50] M. Raidal, V. Vaskonen, and H. Veermäe, J. Cosmol. Astropart. Phys. 09 (2017) 037.

[51] Y. Ali-Haïmoud, E. D. Kovetz, and M. Kamionkowski, arXiv:1709.06576 [Phys. Rev. D (to be published)].

[52] Y. Ali-Haïmoud and M. Kamionkowski, Phys. Rev. D **95**, 043534 (2017).

[53] T. D. Brandt, Astrophys. J. **824**, L31 (2016).

[54] S. M. Koushiappas and A. Loeb, Phys. Rev. Lett. **119**, 041102 (2017).

[55] D. Gaggero, G. Bertone, F. Calore, R. M. T. Connors, M. Lovell, S. Markoff, and E. Storm, Phys. Rev. Lett. **118**, 241101 (2017).

[56] Y. Inoue and A. Kusenko, J. Cosmol. Astropart. Phys. 10 (2017) 034.

[57] M. A. Monroy-Rodríguez and C. Allen, Astrophys. J. **790**, 159 (2014).

[58] B. Carr, M. Raidal, T. Tenkanen, V. Vaskonen, and H. Veermae, Phys. Rev. D **96**, 023514 (2017).

[59] F. Kühnel and K. Freese, Phys. Rev. D **95**, 083508 (2017).

[60] A. M. Green, Phys. Rev. D **94**, 063530 (2016).

[61] K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada, and T. T. Yanagida, Phys. Rev. D **96**, 043504 (2017).

[62] S. Young, C. T. Byrnes, and M. Sasaki, J. Cosmol. Astropart. Phys. 07 (2014) 045.

[63] K. Kohri, T. Nakama, and T. Suyama, Phys. Rev. D **90**, 083514 (2014).

[64] R. Saito and J. Yokoyama, Phys. Rev. Lett. **102**, 161101 (2009); **107**, 069901(E) (2011).
[65] R. Saito and J. Yokoyama, Prog. Theor. Phys. **123**, 867 (2010); **126**, 351(E) (2011).

[66] T. Nakama, J. Silk, and M. Kamionkowski, Phys. Rev. D **95**, 043511 (2017).

[67] T. Nakama, T. Suyama, and J. Yokoyama, Phys. Rev. D **94**, 103522 (2016).

[68] Y. B. Zel'Dovich, JETP Lett. **14**, 180 (1971).

[69] C. W. Misner, Phys. Rev. Lett. **28**, 994 (1972).

[70] A. A. Starobinsky, Zh. Eksp. Teor. Fiz. **64**, 48 (1973) [Sov. Phys. JETP **37**, 28 (1973)].

[71] T. J. M. Zouros and D. M. Eardley, Ann. Phys. (N.Y.) **118**, 139 (1979).

[72] S. L. Detweiler, Phys. Rev. D **22**, 2323 (1980).

[73] A. B. Gaina and I. M. Ternov, Izv. Vuz. Fiz. **31**, 71 (1988) [Sov. Phys. J. **31**, 830 (1988)].

[74] A. Arvanitaki, M. Baryakhtar, and X. Huang, Phys. Rev. D **91**, 084011 (2015).

[75] A. Arvanitaki, M. Baryakhtar, S. Dimopoulos, S. Dubovsky, and R. Lasenby, Phys. Rev. D **95**, 043001 (2017).

[76] K. Kannike, L. Marzola, M. Raidal, and H. Veermäe, J. Cosmol. Astropart. Phys. 09 (2017) 020.

[77] G. Ballesteros and M. Taoso, arXiv:1709.05565 [Phys. Rev. D (to be published)].

[78] F. Takahashi, Phys. Lett. B **727**, 21 (2013).

[79] K. Harigaya, M. Ibe, and T. T. Yanagida, Phys. Rev. D **89**, 055014 (2014).

[80] K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada, and T. T. Yanagida, Phys. Rev. D **95**, 123510 (2017).

[81] D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer, and E. L. Wright, Astrophys. J. **473**, 576 (1996).

[82] K. Inomata, M. Kawasaki, and Y. Tada, Phys. Rev. D **94**, 043527 (2016).

[83] D. Jeong, J. Pradler, J. Chluba, and M. Kamionkowski, Phys. Rev. Lett. **113**, 061301 (2014).

[84] T. Nakama, T. Suyama, and J. Yokoyama, Phys. Rev. Lett. **113**, 061302 (2014).

[85] A. Kogut et al., J. Cosmol. Astropart. Phys. 07 (2011) 025.

[86] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D **81**, 104019 (2010).

[87] A. Barnacka, J. F. Glicenstein, and R. Moderski, Phys. Rev. D **86**, 043001 (2012).

[88] P. W. Graham, S. Rajendran, and J. Varela, Phys. Rev. D **92**, 063007 (2015).

[89] H. Niikura, M. Takada, N. Yasuda, R. H. Lupton, T. Sumi, S. More, A. More, M. Oguri, and M. Chiba, arXiv:1701.02151.

[90] K. Griest, A. M. Cieplak, and M. J. Lehner, Phys. Rev. Lett. **111**, 181302 (2013).

[91] P. Tisserand et al. (EROS-2), Astron. Astrophys. **469**, 387 (2007).

[92] V. Poulin, P. D. Serpico, F. Calore, S. Clesse, and K. Kohri, Phys. Rev. D **96**, 083524 (2017).

[93] L. Lentati et al., Mon. Not. R. Astron. Soc. **453**, 2576 (2015).

[94] R. M. Shannon et al., Science **349**, 1522 (2015).

[95] Z. Arzoumanian et al. (NANOGrav), Astrophys. J. **821**, 13 (2016).

[96] C. J. Moore, R. H. Cole, and C. P. L. Berry, Classical Quantum Gravity **32**, 015014 (2015).

[97] G. Janssen et al., Proc. Sci., AASKA14 (2015) 037, [arXiv:1501.00127].

[98] C. J. Hogan and M. J. Rees, Phys. Lett. B **205**, 228 (1988).

[99] E. W. Kolb and I. I. Tkachev, Phys. Rev. Lett. **71**, 3051 (1993).

[100] E. W. Kolb and I. I. Tkachev, Phys. Rev. D **49**, 5040 (1994).

[101] I. I. Tkachev, Phys. Lett. B **261**, 289 (1991).

[102] E. W. Kolb and I. I. Tkachev, Phys. Rev. D **50**, 769 (1994).

[103] P. Tinyakov, I. Tkachev, and K. Zioutas, J. Cosmol. Astropart. Phys. 01 (2016) 035.

[104] P. Brax, J.-F. Dufaux, and S. Mariadassou, Phys. Rev. D **83**, 103510 (2011).