Axion Models with High Scale Inflation

We revisit the cosmological aspects of axion models. In the high-scale inflation scenario, the Peccei-Quinn (PQ) symmetry is likely to be restored during/after inflation. If the curvature of the PQ scalar potential at the origin is smaller than its vacuum expectation value; for instance in a class of SUSY axion models, thermal inflation happens before the radial component of the PQ scalar (saxion) relaxes into the global minimum of the potential and the decay of saxion coherent oscillation would produce too much axion dark radiation. In this paper, we study how to avoid the overproduction of axion dark radiation with some concrete examples. We show that, by taking account of the finite-temperature dissipation effect appropriately, the overproduction constraint can be relaxed since the PQ scalar can take part in the thermal plasma again even after the PQ phase transition. We also show that it can be further relaxed owing to the late time decay of another heavy CP-odd scalar, if it is present.


Introduction
The axion a appears in the Peccei-Quinn (PQ) mechanism to solve the strong CP problem in QCD [1,2]  where α s is the QCD fine-structure constant, F a is the decay constant of axion, and F a ρσ is the field-strength tensor of gluon. The axion predicts various effects on cosmology [3]. In particular, it is a good candidate of cold dark matter (CDM) in the universe [4,5]. If the PQ symmetry is broken during inflation, the axion obtains isocurvature fluctuations [6][7][8] and the high-scale inflation indicated by the BICEP2 [9] yields too large isocurvature fluctuation inconsistent with observations [10,11]. 1 A more natural situation is that the PQ symmetry is restored during inflation. Actually, for example, the PQ field can have a (positive) Hubble mass through the Planck-suppressed coupling to the inflaton, which easily stabilizes the PQ scalar at the origin during inflation. Then there is no axion isocurvature problem. Instead, the PQ symmetry breaking after inflation produces topological defects: axionic strings and axionic domain walls. The formation of stable domain walls cause cosmological disaster, hence the domain wall number must be equal to one so that the string-wall system collapses soon after the formation [12]. 2 It is shown that the axions radiated from the string-wall system gives a dominant contribution to the relic axion CDM density and it constrains the PQ scale as [14] F a (2.0 − 3.8) × 10 10 GeV. (1.2) If the curvature of the potential of the PQ scalar field at the origin is smaller than its vacuum expectation value (VEV) times its coupling to particles in thermal bath, the thermal inflation is likely to take place before the PQ phase transition. For instance, in supersymmetry (SUSY), the radial component of the PQ scalar is often massless in the SUSY limit and can obtain a soft SUSY breaking mass, which we call saxion. In a class of SUSY axion models, if the PQ symmetry is restored during inflation, the saxion is trapped at the origin due to thermal effects and it causes a brief period of late-time inflation, called thermal inflation [15][16][17][18][19][20], before the PQ phase transition. Even in the non-SUSY case, thermal inflation may take place if the self coupling constant of the PQ scalar is smaller than its coupling to particles in thermal bath. (We call the radial component of the PQ scalar as saxion even in the non-SUSY case.) After thermal inflation ends, the saxion coherent oscillation dominates the universe and finally the saxion decay reheats the universe.
However, this scenario does not always work. This is because the saxion often dominantly decays into the axion pair, resulting in the axion dominated universe, which contradicts with observational constraint on the amount of dark radiation.
In this paper, we study how to avoid the overproduction of axion dark radiation with some examples. We show that, by taking account of the thermal dissipation effect on the saxion [21][22][23], interactions with the thermal plasma can efficiently dissipate the energy of the saxion coherent oscillation into thermal plasma without producing too much axionic dark 1 See, however, Sec. 4 for more on this issue. 2 One can introduce an explicit PQ breaking term to make domain walls unstable even in the case where the domain wall number is larger than one. However, it changes the potential minimum of the axion and the θ angle reappears. In order to make domain walls harmless while the θ angle remains small, some level of tuning is required [13]. We do not consider such a case in this paper. radiation: namely the PQ scalar can participate in the thermal plasma again even after the PQ phase transition. In addition, we also point out that if there exists a heavy CP-odd scalar, whose existence depends on the stabilization mechanism of the saxion potential, the overproduction constraint can be further relaxed owing to its late time decay.
On the other hand, in a class of SUSY axion models, if the PQ scalar participates in the thermal plasma after the PQ phase transition, the axino that might be protected by the (approximate) R-parity will be thermally populated, and it could potentially cause the overproduction problem depending on the stabilization of saxion. We show that the axino overproduction can be avoided with some concrete models.

Decay and dissipation rates in thermal bath
In this paper we consider the hadronic axion model [24]. We discuss the saxion dynamics with the following Lagrangian where Ψ andΨ are PQ quarks with (anti-)fundamental representations of color SU (3), and V stab denotes the potential term that stabilizes the saxion φ at appropriate scale 〈|φ|〉 ≡ f a = F a / 2. We take f a ∼ 10 9−10 GeV to avoid the cosmological and astrophysical constraints. In this section, we concentrate on rather general aspects of the saxion dynamics which does not depend on the detail of V stab . In the next section we will discuss V stab . Let us suppose that φ is placed at the origin during inflation due, for example, to the Hubble-induced mass term. Note that even if φ is placed far from the origin during inflation, it is eventually trapped at the origin φ = 0 due to the efficient particle production at the enhanced symmetry point and its back reaction onto φ [23] unless the coupling, λ, or the reheating temperature, T R , is suppressed. (See also Sec. 4.) Since Ψ andΨ are massless there, they are thermalized after the reheating. It generates the thermal mass for φ as m T ∼ λT . Thermal inflation takes place at the temperature T end < T < T beg where T beg ∼ m φ f a and T end ∼ m φ /λ. 3 After thermal inflation ends, the saxion begins a coherent oscillation around the potential minimum. The perturbative decay rate of the saxion into the axion pair is given by [25] where m s is the saxion mass around the potential minimum, which is same order of m φ . In the following, we do not distinguish m s from m φ unless otherwise stated since it depends on the stabilization of PQ scalar potential. It also decays into the gluon pair, but such a decay mode is subdominant because such a process is loop-suppressed. Therefore the universe would be dominated by the axion dark radiation if the saxion perturbatively decays into the axion pair. Fortunately, the axion overproduction constraint may be relaxed by appropriately taking account of thermal dissipation effect on the saxion in the thermal plasma [21][22][23]. The dominant dissipation comes from the interaction with the thermal plasma via the dimension five operator suppressed by f a , which is obtained from integrating out Ψ andΨ around |φ| = f a .
Then, as one can guess from the dimensional analysis, the reaction rate with the thermal plasma may be proportional to T 3 / f 2 a [26]. In fact, for m φ ≫ g 2 s T , the Landau cut contribution also becomes important, which corresponds to the inverse processes of thermal saxion production [27] proportional to T 3 / f 2 a . 4 On the opposite limit, the pole contribution regulated by its thermal width gives the dominant contribution, and the detailed resummation leads to the following dissipation rate of the saxion coherent oscillation [28,29]: where b is a numerical constant. 5 In the case of quark gluon plasma, b ≃ 1/(32π 2 log α −1 s ) [29]; we adopt this value of b in our numerical calculation. 6 The axions, which are produced non-thermally by the decay of the saxion oscillation, also interact with the thermal plasma. Since their typical energy/momentum are given by m φ that may be much smaller than the cosmic temperature at their production, the dissipation rate of non-thermally produced axions might be different from axion thermal production rate studied in Refs. [30,31]. For p g 2 s T , where p is the typical momentum of the axion, the Landau cut contribution which corresponds to the axion thermal production may be suppressed. Rather, the pole contribution regulated by its thermal width may give the dominant contribution as in the case of the above saxion oscillation. However, contrary to the saxion case, the interaction term can be expressed as aF a µνF aµν = a∂ µ K µ ∼ (∂ µ a)K µ , and hence the axion energy/momentum times temperature should be picked up rather than the thermal mass squared of gauge bosons. As a result, we expect a suppression factor, p 2 /(g 4 s T 2 ), for the nonthermally produced axion dissipation rate compared with the saxion dissipation rate. Using one-loop computation with the Breit-Wigner approximation for the spectral function of gauge bosons, we checked that the suppression factor exists. (See Appendix A.) Our estimation of the axion dissipation rate is given by with a dec and a(t) being the scale factor at the cosmic time Γ −1 φ→2a and at the time t, respectively. In addition, we note here that Eq. (2.4) merely provides an order-of-estimate of the dissipation rate of axion. In order to take account of the uncertainty, we explicitly introduce the parameter C ∼ (1).

Boltzmann equations
After the thermal inflation, the energy densities of the saxion oscillation, axion, and radiation, which are denoted as ρ φ , ρ a , and ρ r , respectively, evolve aṡ where H is the expansion rate of the universe. Thus, after a few Hubble times after the PQ phase transition, due to the dissipation effect, the temperature increases to 18 GeV being the reduced Planck scale. Notice that T max ∼ T c when the dissipation rate is smaller than the expansion rate for T ∼ T c so that the effect of the dissipation is ineffective. If the dissipation rate can become larger than the Hubble parameter, the energy density of the saxion coherent oscillation is expected to be transferred to the radiation within a few Hubble time. In other words, the PQ scalar takes part in the thermal plasma again. The condition is written as In the above estimation, we have implicitly assumed that the saxion dominates the universe before the dissipation becomes effective. However, in general, it is possible that the saxion dominantly decays into axions at first, and then the axions which dominate the universe transport their energy into radiation. The decay mode into axion dominates when Note that the ratio between the saxion decay rate into the axion pair and the Hubble parameter H at the end of thermal inflation is given by (2.13) If this ratio is larger than one, and also if Γ φ→2a is larger than Γ (dis) φ , the saxion soon decays into the axion pair within one Hubble time. In this case, the temperature can be estimated as where δ is a numerical factor; The condition for axion thermalization is given by If these conditions [Eqs. (2.11), (2.15)] are satisfied, the temperature increases up to T ∼ T max within a few Hubble time after thermal inflation and the saxion coherent oscillation disappears without producing too much axion dark radiation. As one can see, in the case of the axion CDM, the mass scale m φ should be larger than (1) PeV in order for the PQ scalar to take part in the thermal plasma, 7 which implies the maximum temperature to be T max ∼ (10 7 ) GeV for m φ ∼ (1) PeV. Note that such a value of m φ is suggested by high-scale SUSY breaking models [32][33][34]. The detailed discussion on which parameters the overproduction of axion dark radiation can be avoided is performed with some examples in the following, since it is model dependent.
Before going into details of models, let us briefly give general comments. PQ quarks, Ψ and Ψ, become massive after the saxion obtains a VEV and they can soon decay if there is a mixing between PQ quarks and SM fermions [23]. In a class of SUSY axion models, for (m φ , f a ) ∼ ( (1) PeV, (10 10 ) GeV), gravitinos are thermally produced [35], but its abundance is not so large and it does not cause cosmological problems for m 3/2 ∼ (100) TeV. The axion multiplet, saxions, axions and axinos, will also be thermally populated if the dissipative effect is efficient. In particular, the axino whose stability might be protected by the (approximate) R-parity must decay well before Big-Bang nucleosynthesis (BBN) to avoid cosmological constraints. The axino mass crucially depends on the saxion stabilization model. Below we will discuss it with two models (Model 2 and 3).

Cosmology: model-dependent discussion
In the following subsections, we study in detail the dynamics of the saxion after the thermal inflation, and show how the overproduction constraint of axion dark radiation is relaxed with some bench mark models. First, we consider the simple non-SUSY case (Model 1), where it is shown that the overproduction constraint is relaxed by interactions with the thermal plasma. Then, we consider SUSY axion models with two stabilization mechanisms of the saxion potential: radiative stabilization (Model 2) and higher dimensional superpotential with another PQ field (Model 3). In particular, in the latter case (Model 3), we show that the overproduction constraint is further relaxed by the late time decay of the heavy CP-odd scalar.
In addition, in a class of SUSY axion models, there exists the axino which might be protected by the (approximate) R-parity, and hence it could potentially overclose the universe due to the thermal production, if the dissipation of the PQ scalar is efficient. Since the axino overproduction problem crucially depends on the stabilization mechanism of the saxion, we also discuss how to avoid this problem separately with two cases (Model 2 and 3).

Model 1
Let us first consider the non-SUSY case where the saxion potential is simply given by (3.1) In this case, all the necessary ingredients are already given in the previous subsection. Hence one can calculate the effective neutrino number of axion dark radiation for each parameter (m φ , f a ) by numerically solving Eqs. (2.6) -(2.8).
; with ρ • being the energy density of •. Note that ρ φ and ρ a denote the energy density of the saxion coherent oscillation and that of the axion produced from it respectively, and both does not contain that from thermal plasma. See also footnote 7.
The results are shown in Fig. 1, which shows evolution of various quantities after the PQ phase transition as a function of Hubble time H −1 normalized by that at the PQ phase transition, H −1 PT . The black, pink and blue lines represent evolution of the energy densities of the saxion coherent oscillation, radiation, and the non-thermal axion produced from the saxion coherent oscillation normalized by that of the initial saxion coherent oscillation, ρ φ,PT . In the case with (m φ , f a ) = (50 PeV, 10 10 GeV) [Left Panel of Fig. 1(a)], the saxion coherent oscillation first loses its energy into axion and radiation within a few Hubble times. Then, the non-thermally produced axions interact with radiation and reduce their number. In this specific parameter, the predicted effective number of extra radiation is ∆N eff ∼ 0.4. On the other hand, in the case with (m φ , f a ) = (10 TeV, 5 × 10 8 GeV) [Right Panel of Fig. 1(b)], the saxion coherent oscillation can lose almost all the energy before the decay into axions dominates. The resultant axion dark radiation is ∆N eff ∼ 0.2. Fig. 2 shows contour of ∆N eff = 1 on (m φ , f a ) plane for T end = m φ (corresponding to λ ∼ 1) (left panel) and T end = 10m φ (corresponding to λ ∼ 0.1) (right panel). 8 In the figure, we take account of the fact that our estimation of the dissipation has uncertainty factor owing to both the model dependence (e.g., contributions from other gauge group, gauginos) and theoretical uncertainties. Thus, we vary the dissipation rate and show how much the bound changes; for this purpose, we have varied the C-parameter in Eq. (2.4). The upper and lower boundaries of the band shown in Fig. 2 correspond to C = 1 and C = 10, respectively. Above the band, ∆N eff becomes smaller than one. The black dashed line represents the contour ∆N eff = 1 without the thermal dissipation. We can see that the region with small enough ∆N eff is significantly enlarged by taking account of the effect of axion dissipation. In particular, in the case where φ couples to the PQ fermions relatively strongly (i.e, λ ∼ 1), m φ is required to be close to f a if the effect of thermal dissipation is neglected; this is in order to avoid the saxion domination after the PQ phase transition. With the proper inclusion of the effects of thermal dissipation, m φ much smaller than f a becomes allowed even if λ ∼ 1.

Model 2
Next, let us consider the SUSY axion model. First, we focus on the possibility that the running of the soft mass of φ induces the radiative symmetry breaking at an appropriate scale [36,37]. We consider the following superpotential: where φ denotes the PQ scalar,Q(Q) and L(L) are chiral multiplets in the (anti-)fundamental representation of SU (5). Then, the renormalization group equations (RGEs) for the Yukawa coupling constants are where t = log E (with E being the energy scale) and (3.4) In addition, the renormalization group equations for the soft SUSY breaking parameters are given by where m X (with X = Q,Q, L,L) denote the soft SUSY breaking scalar mass parameters for the scalar components in corresponding chiral multiplet. Here we have ignored gaugino masses by assuming that gaugino masses are much smaller than those of sfermions as in the pure gravity mediation model. In this model, the dynamics of the PQ scalar after the PQ phase transition is almost the same as the model 1. Therefore, from Fig. 2, we can understand the parameter region on (m φ , f a ) plane, where the overproduction of axion dark radiation is avoided. For instance, the axion CDM implies that the mass scale, m φ , should be larger than (1) PeV. To see if the radiative symmetry breaking really occurs, we have solved the RGEs with the boundary condition at the GUT scale as m 2 φ = m 2 Q = m 2 L = (10 PeV) 2 and λ Q = λ L = 1.0 and 0.7. The result is shown in Fig. 3. It is seen that m 2 φ becomes negative at some intermediate scale. In particular, for λ Q,L ≃ 0.7, the VEV of φ can be around 10 10 GeV. Thus the cosmological scenario studied in the previous section works. However, this model suffers from the axino overproduction. In this model, the axino mass arises only radiatively and is suppressed compared with the gravitino mass [37,38]. It is estimated as mã ∼ (5/16π 2 ) 2 λ 4 m 3/2 . Thus it is lighter than the gauginos and it likely becomes the lightest SUSY particle (LSP). If the thermal dissipation is effective, the axion multiplet takes part in the thermal plasma and hence axinos are thermally populated, which results in overproduction of the axino. Hence the presence of stable axino LSP is problematic. Below we list some possibilities to avoid the axino overproduction.
• If one introduces another LSP, the problem may be avoided. For example in a singlet extensions of the minimal SUSY standard model [39], the singlino-like neutralino can be the LSP. In this case the axino can decay into the singlino LSP through the singlino-gaugino mixing, which is suppressed by the SUSY breaking scale. The axino decay temperature is then estimated as where M denotes gaugino mass scale and µ is the higgsino mass. In order for the axino to decay well before BBN, relatively light neutralino mass spectrum is needed. In addition, singlino also has to decay before BBN, which requires (small) R-parity violation. Note that if the axino decays after the QCD phase transition, the axion abundance is diluted and upper bound on the PQ scale is relaxed [40].
• The axino may become heavier if there is a sizable A-term, which yields the axino mass as mã ∼ (λ 2 /16π 2 )A φ [38] and it can exceed the anomaly mediation contribution to the gaugino masses [41] if A φ ∼ m 3/2 . Actually we naturally have A φ ∼ m 3/2 if the SUSY breaking field z has no gauge quantum number, as in the Polonyi model. The gauginos are also expected to obtain masses of ∼ m 3/2 , but the axino can be heavier if the coupling between gaugino and Polonyi field happens to be relatively small. Then the axino decays into gaugino. Note that, in this case, the assumption on RGEs, Eqs. (3.3)-(3.5), is marginally satisfied since the gauginos are lighter than the gravitino. If the Wino is LSP, Winos produced by the axino decay annihilate efficiently and its abundance can become smaller than the DM abundance [42]. In this case, however, there arises a cosmological Polonyi problem [43,44] which cannot be solved even in the presence of thermal inflation. In Sec. 3.4, we will discuss this issue and how to solve it.

Model 3
Let us consider the following model with two PQ scalar φ andφ [45]: where φ andφ are assumed to have PQ charges of +1 and −n, respectively, and M is a cutoff scale. Then the scalar potential is (3.9) where the parameters m 2 φ and m 2φ are both taken to be positive. For computational simplicity, we assume m 2φ ≪ m 2 φ . Then the potential minimum lies at In this model, the fermionic component of φ andφ, which we denote byã andã and call them as axino, obtain both Dirac and Majorana masses after PQ scalars get VEVs. Then the axino masses are generated as mã ∼ m φ . Since they are as heavy as the gravitino, they decay into the gluino and gluon: (3.11) The axino decay temperature is then estimated as Thus, they are harmless. Note that the massless axion is a mixture of the angular component of φ andφ. The real components of φ andφ are also mixed with each other and both φ andφ begins a coherent oscillation after thermal inflation ends [46]. Although the dissipation effect acts only on φ, the time scale of the mixing (∼ m −1 φ ) is much faster than that of the cosmic expansion and hence both the energy density of φ andφ are efficiently dissipated.
Importantly, there also exists a heavy CP-odd scalar a ′ with a mass ∼ m φ as a mixture of the angular component of φ andφ. After the end of thermal inflation, the heavy CP-odd scalar may also begin an oscillation with an amplitude ∼ f a . Its dissipation rate is expected to be the same as that of the axion except for the redshift factor, where x = m φ /(g 4 s T ). Hence above the contour shown in Fig. 2 it is expected to be dissipated away into the thermal plasma. On the other hand, below this contour the heavy CP-odd scalar survives from interactions with the thermal plasma, and later it decays into gluons/gluinos through one-loop processes. Since the entropy production becomes more significant for a smaller mass, m φ , the overproduction of the axion dark radiation can be avoided in all the region shown in Fig. 2 owing to the combination of the thermal dissipation and of the late time entropy production from the heavy CP-odd scalar. Fig. 4 shows evolution of same quantities in Fig. 1 plus the energy density of heavy CP-odd scalar after the PQ phase transition with (m φ , f a ) = (50 PeV, 10 10 GeV) and (1 TeV, 10 10 GeV) in the left and right panel respectively. Here we have simply assumed that the initial energy density of heavy CP-odd scalar is the same as that of saxion. As one can see from the right panel [ Fig. 4(b)], the late time decay of the heavy CP-odd scalar successfully dilutes the axion dark radiation.

Dilution of unwanted relics and baryon asymmetry
Thermal inflation dilutes the preexisting unwanted relics such as the moduli, Polonyi and gravitinos as well as the baryon number. If we consider the case where the thermal dissipation is efficient, f a is relatively small while m φ is of the order of PeV or larger. Then, the e-folding number of the thermal inflation is not so large. In the present case, the dilution factor ∆ is estimated to be , with ρ • being the energy density of •. Note that ρ φ and ρ a denote the energy density of the saxion coherent oscillation and that of the axion produced from it respectively, and both does not contain that from thermal plasma. See also footnote 7.
For example, the final baryon asymmetry after thermal inflation is evaluated by multiplying the original baryon asymmetry by the inverse of this factor: where (n B /s) reh denotes the baryon-to-entropy ratio evaluated as if there were no thermal inflation. Since the dilution factor is not so huge in the typical parameters for which the dissipative effect is efficient, it is possible that the preexisting baryon asymmetry, created e.g. by the Affleck-Dine mechanism [47,48], survives the dilution due to thermal inflation and it explains the present baryon asymmetry of the universe. Note also that the temperature after thermal inflation increases to T max ∼ 10 7−8 GeV, and hence it is also possible that thermal leptogenesis [49] works if there is a mild degeneracy among right-handed neutrino masses [50,51]. On the other hand, the above dilution factor is not sufficient to solve the cosmological Polonyi/moduli problem. If there is a singlet Polonyi/moduli field, its abundance is given by where T (inf) R is the reheating temperature after inflation and z i is the initial amplitude of the Polonyi field. If the Polonyi mass is of the order of the gravitino, its decay produces too much LSPs. Actually, in the model 2 discussed above, we may need a singlet Polonyi field in order to make the axino heavy. In this case the cosmological Polonyi problem is not solved by thermal inflation, and we need to introduce a (small) R-parity violation to make the LSP unstable or to rely on the adiabatic suppression mechanism [52,53] to suppress the Polonyi abundance. In the model 3, on the other hand, we do not need such a singlet Polonyi field and hence there is no Polonyi problem. 9

Discussion
In this paper we have revisited the cosmology of axion models. If the curvature of the potential of the PQ scalar field is smaller than the PQ scale, like in a class of SUSY axion models, it is likely that the PQ scalar field is trapped at the origin and it causes thermal inflation. We have discussed several possibilities to avoid the overproduction of the axion dark radiation from the saxion decay after the thermal inflation. In general, the overproduction constraint is relaxed by taking thermal dissipative effects on the saxion coherent oscillation into account. This is because, owing to the interaction with the thermal plasma, the PQ scalar can participate in the thermal plasma even after the PQ phase transition. This scenario works for the mass scale, m φ , larger than (10) TeV -(10) PeV for f a ∼ 10 9 GeV -10 10 GeV. We have demonstrated it with some explicit models of the saxion stabilization. In particular, we have shown that if there exists a heavy CP-odd scalar, depending on the stabilization model, it can dilute the axion dark radiation due to its late time decay. As a result, the overproduction constraint is further relaxed. Those studies have significant implications to the case of high-scale inflation. However, we note here that our results are not restricted to such a case.
One interesting prediction of the SUSY axion model is that it naturally deforms the gaugino masses [36,56]. If the gaugino masses are dominated by the anomaly mediation effect, the additional contribution from the PQ multiplet is generally the same order. Therefore, if gauginos are found and their masses are measured at the LHC, it may indicate a signal of the axion model.
Hereafter we comment on the recent claim of the discovery of the B-mode polarization in the cosmic microwave background anisotropy by the BICEP2 experiment [9], which indicates the high-scale inflation, such as the chaotic inflation [57]. The derived inflation scale in terms of the Hubble scale is as high as H inf ≃ 10 14 GeV.
A prediction of the present scenario is the suppression of inflationary gravitational waves (GWs) at the high frequency range at which future space laser interferometers have high sensitivities [58]. The energy density of inflationary GWs with frequency f > f c where is suppressed by the factor ∆ 4/3 and hence it may be difficult to detect [59] if the dilution factor ∆ is sizable. Instead, if the PQ phase transition is first order, there is a possibility that GWs from bubble collisions will be in the observable range [60]. So far we have focused on the hadronic axion model, because it is the simplest setup for giving the domain wall number one. In the DFSZ axion model [61,62], the PQ scalar couples to the Higgs multiplets. However, in this model the domain wall number is larger than one, and hence there is a serious domain wall problem if the PQ symmetry is restored during inflation. In order to make the domain wall number one, we need to introduce five pairs of additional PQ quarks. Thus the coupling to the PQ quarks leads to thermal inflation as studied in this paper. An important feature in this model is that the saxion can dominantly decay into the Higgs boson pair. Hence we do not need thermal dissipation effects to avoid the axion overproduction.
Finally we comment on the case in which the PQ symmetry is broken during inflation. In such a case, the high scale inflation may be excluded because of too large axion isocurvature perturbation. This argument has some loopholes. During inflation, the saxion may have a field value much larger than f a if the saxion obtains the negative Hubble induced mass term. Then the magnitude of axion isocurvature fluctuation is highly suppressed [63]. The fate of the saxion field with such a large initial amplitude strongly depends on the saxion stabilization model. In the model 2, the saxion can have a field value of ∼ M P during inflation, but it is eventually trapped at the origin after the onset of coherent oscillation as extensively studied in Ref. [23]. 10 The subsequent dynamics is the same as that studied in this paper. In the model 3, the saxion tracks the temporal minimum |φ| ∼ (H M n−2 ) 1/(n−1) as the Hubble parameter decreases and it relaxes to the true minimum without restoration of the PQ symmetry if the coupling, λ, or the reheating temperature, T (inf) R , is suppressed [46]. For n = 3, for example, the magnitude of CDM isocurvature perturbation is given by , (4.2) where θ i is the initial misalignment angle and ρ a and ρ CDM are the energy density of axion and CDM, respectively. 11 Since the Planck data combined with the WMAP9 polarization data gives an upper bound as S CDM 1×10 −5 [66], the isocurvature constraint can be avoided. However, the axion cannot be the dominant component of CDM in this case. See also Refs. [67][68][69] for recent discussion in this direction. The saxion coherent oscillation with an amplitude of ∼ f a remains in this case, but it is not cosmologically problematic [46].

A One-loop estimation of soft axion dissipation
In this appendix, we estimate the thermal dissipation rate of soft relativistic axions with p ≪ g 2 s T at the one-loop order. In principle the following estimation may be made more precise by the resummation of many diagrams, but we postpone this task as a future work.
The dissipation rate of relativistic axion is given by where q 0 = p 0 − k 0 , d a = N 2 c − 1 = 8, f B is the Bose-Einstein distribution function, and ρ T /L denotes the spectral function for the transverse/longitudinal mode of gauge boson in the thermal plasma. Roughly speaking, for p 0 = p ≫ g 2 s T , the cut contribution of spectral function becomes important since the typical center of mass energy for scattering processes becomes larger than thermal mass, pT ≫ g 2 s T 2 . And for p 0 = p ≫ g s T , the dominant contribution is calculated by means of the Hard Thermal Loop approximation in Ref. [30,31]. On the other hand, for p ≪ g 2 s T , the near pole contribution of spectral function becomes important, which is smeared by interactions in the thermal plasma.
To estimate this integral, we approximate the spectral density for the transverse mode with the Breit-Wigner form: where m ∞ is the asymptotic mass: for instance in the case of non-abelian gauge theory with N F flavor, it is given by m 2 ∞ = [C(ad) + N F T (F)]g 2 s T 2 /6 with C(ad) being the adjoint Casimir and T(F) being the normalization index. Here we estimate the width by large angle scatterings, Γ p ∼ g 4 s T 3 /p 2 [28]. Since the loop-integral is typically dominated by a hard momentum ∼ T , we omit contributions from longitudinal mode ρ L , which is suppressed by the residue [70]. After some calculations, one can express the dissipation rate with The asymptotic behavior of f (x) is f (x) ∼ const. for x ≪ 1 and f (x) ∝ x −2 for x ≫ 1.
Notice that our estimation of Γ (dis) a is expected to be an order-of-estimate of the dissipation rate because it is based on a simple one-loop calculation. Thus, in order to take account of the uncertainty in our result, we introduce the numerical coefficient C and vary it in our numerical analysis.