Inflationary Gravitational Waves and the Evolution of the Early Universe

We study the effects of various phenomena which may have happened in the early universe on the spectrum of inflationary gravitational waves. The phenomena include phase transitions, entropy productions from non-relativistic matter, the production of dark radiation, and decoupling of dark matter/radiation from thermal bath. These events can create several characteristic signatures in the inflationary gravitational wave spectrum, which may be direct probes of the history of the early universe and the nature of high-energy physics.


Introduction
The early universe is a good laboratory for high-energy physics because the temperature can be much higher than the reach of the accelerator experiments. Thus we may be able to probe high-energy physics by observing some relics of the hot early universe. The anisotropy of the cosmic microwave background (CMB) carries rich information on the primordial density perturbation, which is considered to be generated during inflationary era. Hence the precise observation of CMB gives a clue to inflation models; it is a great success of cosmology for probing high-energy physics.
However, it is rather unknown what have happened in the era between the inflation and the big-bang nucleosynthesis (BBN). To go beyond, one of the promising such relics that carry direct information on the early universe is the gravitational waves (GWs), since GWs propagate without interfered by matter and radiation. In this paper, we focus on inflationary GWs as a possible probe of the early universe.
In the presence of decaying matter into dark radiation, the GW spectrum is deformed in a nontrivial way [28].
In a realistic setup, the GW spectrum would be affected in a complicated way. For example, in the Peccei-Quinn (PQ) model for solving the strong CP problem [29][30][31], the PQ phase transition may have occurred at around the cosmic temperature T ∼ 10 9 -10 12 GeV, and huge amount of entropy could have been released in association with the phase transition. Relativistic axions may have been also efficiently produced by the decay of the PQ scalar condensate, which would contribute as dark radiation. All these phenomena significantly affect the spectrum of inflationary GWs. In other words, it may be possible to access the high-energy phenomena and underlying physical processes by studying the detailed spectrum of inflationary GWs.
In this paper we study the inflationary GW spectrum in detail in various setups having concrete models in mind. The combinations of above mentioned effects make several characteristic features in the spectra, which enable us to infer information on the early universe phenomena.
First we review basic properties of GWs in section 2, including the GW spectrum and effects of the anisotropic stress. In section 3, we exhibit some simple examples of the GW spectrum as a first step to the following section. In section 4, we combine all these effects to show the realistic GW spectra with several typical examples. In section 5 particle physics motivated models will be provided in which the GW spectrum actually becomes complicated. Section 6 is devoted to conclusions.

Background evolution
Before discussing the evolution of GWs we first explain that of background, which we assume throughout this paper to be the Friedmann-Robertson-Walker (FRW) universe with negligible curvature. The FRW metric and its tensor perturbation are given by where a(t) is the scale factor and h ij = h ji satisfies transverse and traceless condition: The evolution of the FRW universe is described by the Friedmann equation where H =ȧ/a, and M P ≃ 2.4 × 10 18 GeV is the reduced Planck mass. The total energy density of the universe, ρ tot , generally include where ρ vac , ρ m , ρ r and ρ X represent the vacuum, matter, (visible) radiation and dark radiation energy densities, respectively. 1 The radiation energy density ρ r is related to the cosmic temperature T as

4)
where g * is the relativistic degrees of freedom. Unless otherwise stated, we use g * = 228.75, which is the value in the minimal supersymmetric (SUSY) standard model (MSSM) at high enough temperature. 2 The amounts of other components depend on the particle physics model and there can be energy transfers among these components. Thus the Hubble expansion rate reflects their behavior in the early universe and it directly affects the GW spectrum as shown below. We will study concrete setups in the following sections.
2.2 Evolution of GWs: case without dark radiation

Evolution equation
Substituting the metric (2.1) into the Einstein equation, we get the equation of motion of GWs:ḧ where Π ij is the anisotropic stress of the energy-momentum tensor, which satisfies Π ii = Π ij,i = 0. We decompose h ij using polarization tensors e +,× ij , 1 In this paper, we call the energy density with w = −1 (with w being the equation-of-state parameter) "vacuum energy".
2 Including other fields than in MSSM causes deviation from this value, but we neglect it as small.
Here k = |k|, and we have performed the same decomposition and transformation on Π ij as is done on h ij . Also, we define the polarization tensors so that they satisfy e λ ij = e λ ji , e λ ii = e λ ij,i = 0 and e λ ij e λ ′ * ij = δ λλ ′ . For later use, we introduce the variable u defined by , (2.8) and rewrite (2.7) to get where the prime denotes the derivative with respect to u and H u ≡ a ′ /a. If a(t) ∝ t p as in the radiation dominated (RD) era (p = 1/2) or the matter-dominated (MD) era (p = 2/3), we obtain u = p 1−p k aH . Imposing the initial condition of h(t, k, λ) → h prim (k, λ) for t → 0, the solution of (2.9), when the anisotropic stress is neglected, is given by and where j i are the i-th spherical Bessel function: From this solution, it is easily seen that h(u, k, λ) ∼ const. for the modes outside the horizon (k ≪ aH) and h(u, k, λ) ∝ a −1 for the modes inside the horizon (k ≫ aH).

GW spectrum: modes entering the horizon at the RD era
Here we define basic quantities used in the following sections. The energy density of the GWs is given by (see appendix A for details) where ρ GW (t, k) denotes the energy density of tensor perturbation per logarithmic frequency. In addition, · · · osc denotes the oscillation average (and hence the above expression is relevant only for the sub-horizon mode). Taking the ensemble average, we obtain where P h is the power spectrum of GWs. The GW spectrum Ω GW (t, k) is defined as

JCAP01(2014)040
For the GW modes entering the horizon at the RD era, the present value of Ω GW is evaluated as where r is the tensor to scalar ratio, T hi is the temperature at which the mode enters the horizon and T eq is the temperature at the matter-radiation equality. 3 (Here, we assume that there is no entropy production after the horizon entry.) Ω GW is weakly dependent on k through T hi (k) and the tensor spectral index n t . The relation between the present frequency of GWs and the temperature of the universe at which the corresponding mode entered the horizon is given by (2.17) Since future space-based GW detectors are most sensitive to the GW with frequency around 0.1-1 Hz, we can probe the early universe physics with very high temperature through the GW observations.

GW spectrum: modes entering the horizon at the non-RD era
In the above, we have assumed the RD universe at the horizon entry of GWs. This is not guaranteed in general in the early universe where various forms of fluid can dominate the energy density. Let us suppose that the universe has the equation of state w(> −1/3), the ratio of the energy density to the pressure, at the horizon entry of GWs. Then the GW spectrum scales as In the MD universe (w = 0), Ω GW (t 0 , k) ∝ k −2 .
If there is a short period of inflation, the expression is a bit complicated. Let us suppose that the equation of state changes as w 1 → w 2 → w 3 with w 1 , w 3 > −1/3 and w 2 < −1/3. The calculation is done in appendix B and the result is for the modes entering the horizon at the intermediate regime. If w 1 = w 3 = 1/3 and w 2 = −1, we obtain Ω GW (t 0 , k) ∝ k −4 . If w 1 = 1/3, w 3 = 0 and w 2 = −1, we obtain Ω GW (t 0 , k) ∝ k −6 .

Evolution equation
In the presence of relativistic particles with weak or no interaction, the r.h.s. of (2.9) does not vanish and affects the evolution of GWs [24][25][26][27][28]. Let us consider "dark radiation" X, noninteracting relativistic particles, contributing to the anisotropic stress. 4 The r.h.s. of (2.9) 3 Since at least two species of neutrinos have masses and non-relativistic at present, we use g * at T = Teq rather than that at T = T0 (present temperature) to avoid confusion. At T = Teq, all neutrinos are relativistic. See also appendix A. 4 We call non-interacting relativistic degrees of freedom at the time of our interest as "dark radiation." Thus, the dark radiation in our discussion may be massive, and may not correspond to the dark radiation in the present universe.
-5 -is written as an integration including the metric perturbation h(u, k, λ) (see ref. [28] for derivation), where ρ X is the energy density of X, and j 2 (u) = (3 − u 2 ) sin(u) − 3u cos(u) /u 3 . Also, we have omitted trivial indices. The right-hand side of (2.20) is due to the backreaction of dark radiation on GWs with kernel j 2 (u)/u 2 . Note that j 2 (u)/u 2 is suppressed when u ≫ 1 while h ′ (u) ≃ 0 at u ≪ 1, hence the anisotropic stress affects the evolution of GWs only around u ∼ 1, i.e., around the horizon entry in the RD or MD universe.

Overall normalization
The presence of non-interacting particle X at the time of horizon entry or the last-scattering causes a change in the overall normalization of the GW spectrum as shown in ref. [28]. Here we briefly repeat the result of ref. [28]. First, the amount of dark radiation is written in terms of the effective neutrino species N eff as where ρ ν and ρ γ are the energy densities of the neutrino and photon, respectively. With a standard matter content one expects N (SM) eff = 3.046, and we define The dependence of Ω GW (t 0 , k) on ∆N eff comes from g * and g * s in (2. 16), where, here and hereafter, the superscript "(std)" is for quantities in the case without dark radiation (i.e., ∆N eff = 0). Also, the relation between g * (T eq ) and g Then we define the overall factor C 1 as which is plotted in figure 1.

JCAP01(2014)040
This C 1 factor is partially compensated by the effect of anisotropic stress caused by X. If X exists before the horizon entry, the r.h.s. of (2.20) can be simplified as [24] r.h.s. of (2.20) where f X = ρ X /ρ tot is the energy fraction of X. In this case, assuming the RD universe, the suppression of GW spectrum caused by X is calculated analytically. The suppression factor is defined as 5 whose dependence on f X was analytically derived in refs. [4,25] and plotted in figure 2.
The relation between f X and ∆N eff is given by Then the present GW spectrum reads for the mode entering the horizon at the RD era. This gives the overall normalization of the GW spectrum in the low-frequency limit, k < k EW , where k EW is the comoving Hubble scale at around the electroweak phase transition. The total modification on the overall normalization on Ω GW , C 1 × C 3 , is plotted in figure 1.
In the following sections, we calculate the GW spectrum by using numerical calculation, taking account of the effects of anisotropic stress. However, the normalization related to the C 1 factor is not included in the following calculations because it is model-dependent; C 1 depends whether X remains non-interacting and relativistic until present or not. Thus one should note that Ω GW given in the figures in the following sections should be multiplied by C 1 if the X particle behaves as dark radiation until today.

Illustration with simple examples
As we have mentioned, the spectrum of the inflationary GWs is sensitive to the thermal history of the universe. To see basic features of GW spectrum, in this section, we exhibit some simple examples for the GW spectrum modified by the phase transition, entropy production and dark radiation.

Phase transition
As a first example, let us consider a scalar field φ with a symmetry breaking potential, e.g., We consider the case where φ is trapped at the origin due to thermal effects [23]. If the interaction of φ with the particles in thermal bath is strong enough, the vacuum energy dominates the universe, and a brief period of inflation occurs, as in the case 5 The notations C1 and C3 follow those of ref. [4].  of thermal inflation [19][20][21][22]. The phase transition happens after inflation and the subsequent oscillation of the field is assumed to instantly decay into radiation. The background equations we solve are (2.2) with where t PT is the cosmic time at the phase transition. We numerically solved these with (2.9), and the resultant Ω GW is shown in figure 3. We varied the ratio of ρ r to ρ tot ≡ ρ r + ρ vac at the time of phase transition. The horizontal axis is normalized with k PT , which satisfies k = aH at the phase transition. Note that the ratio of the spectrum in k ≫ k PT to that of k ≪ k PT is equal to ρ r /ρ tot just before the phase transition, which is from the fact that ρ GW (k ≫ aH) ∝ a −4 . GWs inside the horizon at the phase transition are diluted by the newly-produced radiation, while those outside the horizon remains h ij = const. Note the 1 . GW spectrum with phase transition and instant decay into radiation. We assumed that the universe is radiation dominated before the vacuum energy dominates it, and varied the ratio of radiation energy density to the total energy density at the phase transition. characteristic oscillatory feature around k ≃ k PT . The reason for the oscillatory feature and the oscillation period is explained in appendix B. Finally, note that in the limit of long duration of thermal inflation, the spectrum scales as ∝ k −4 as shown also in appendix B.

Entropy production
Next we consider the case of late-time entropy production, i.e., the case where some matter dominates universe, and then it decays into radiation. The background equations we solve are (2.2) andρ where Γ is the decay rate of the non-relativistic matter. The result is shown in figure 4. We varied the ratio of ρ r to ρ tot ≡ ρ r + ρ m at the decay. The horizontal axis is normalized with k decay , which is defined as the comoving Hubble scale aH at t = Γ −1 . Note that the ratio of the spectrum in k ≫ k decay to that of k ≪ k decay is equal to ρ r /ρ tot at the decay for the same reason written in the previous subsection. In this case the spectrum scales as ∝ k −2 for the mode entering the horizon at the non-relativistic matter dominated era. Also note that there is no oscillatory feature in the GW spectrum in contrast to the previous case.

Phase transition and dark radiation
Let us consider the case where vacuum energy of a scalar field dominates the universe as in the case of section 3.1, but at a certain time the energy is instantly converted to dark radiation X. We numerically solved the Friedmann equation (2.2) with 6 . GW spectrum with entropy injection. We varied the ratio of radiation energy density to the total energy density at t = t dec ≡ Γ −1 .
Then we expect the following features in the GW spectrum: • For large wavenumber k ≫ k PT , there is the same suppression as figure 3.
• For small wavenumber k ≪ k PT , there is a suppression due to the anisotropic stress, caused by the r.h.s. of (2.20).
The reason for no suppression by the anisotropic stress at k ≫ k PT is that, for GWs of such large wavenumber, there do not exist X particles at the time of their horizon entry. The results of numerical calculations are shown in figure 5. In the figure we varied Λ so that ∆N eff (after the phase transition) becomes 1, 2, 5 and 100. Note that ∆N eff here is evaluated assuming that X is relativistic and survives until today. If X decays into radiation at some epoch, or if there exists another entropy production, the X abundance at the epochs of the BBN and radiation-matter equality can be reduced, hence ∆N eff ≫ 1 does not necessarily conflict with observations. In such a case the overall normalization of the GW spectrum changes, but the shape of the spectrum does not change. One sees a dip around k ≃ k PT , which could be a smoking-gun signal of the phase transition followed by the production of dark radiation.

Entropy production and dark radiation
Next let us consider the case where some massive particle dominates the universe as in the case of section 3.2, but then it decays into dark radiation X [28]. We numerically solved (2.2) anḋ The results of numerical calculations are shown in figure 6. One also sees a characteristic dip around k ≃ k PT as in the previous case.

Decay of dark radiation into visible radiation
Another interesting possibility is that X is produced at some time or has existed from the beginning, and then decays into visible radiation while X is still relativistic. In this case the background equations we should solve are where m is the (small) mass of X, Γ is the decay rate of X and F X (t, E) is defined so that X with energy E -E + dE carries energy density of F X (t, E)dE at the time t. 6 It satisfies dEF X (t, E) = ρ X (t). We solved (2.2), (2.20), (3.11), and (3.12), varying the initial ratio of ρ X to ρ r and the energy dependence of F X . The result is shown in figure 7 for the case where 6 Without decay, the energy density carried by X with energy E -E + dE in a comoving volume falls proportional to a −1 : the dark radiation initially dominates the universe. We consider three cases with different distribution functions: the line (i.e., F X (t, E) ∝ δ(E 0 )), bosonic, and fermionic thermal distributions. Note that the spectrum has a hill around k ≃ k dec , where k dec for the line distribution is defined as the comoving Hubble scale k = aH at t = (mΓ/E) −1 . For thermal distribution it is defined as k = aH at t = (mΓ/Ẽ) −1 withẼ ≃ T . We found no significant differences in the GW spectrum among these distributions. Thus, in studying the relativistic decay of dark radiation in the next section, we approximate that it has the line spectrum.
If we include decay, the energy density decreases with decay rate Γ suppressed by the γ-factor: Again using E(t), dE(t) ∝ a −1 , we get (3.11).

Examples of the GW spectrum
In this section we study the spectrum of the GW in cases with the combinations of events discussed in the previous section, which may be realized in some models motivated by particle physics. (See the flow chart given in figure 8.) In our study, we perform the analysis as general as possible, without specifying underlying models. Examples of particle-physics models realizing the scenarios in this section will be discussed in the next section. In this subsection we assume that the universe has undergone the following events in a time ordering: 1. A brief period of thermal inflation is caused by a scalar field φ.
2. After the phase transition, φ instantaneously decays into radiation with short mean free path.
3. X particles decouple from the thermal bath, after which X particles behave as dark radiation.
4. Part of the decoupled particles X decays into visible radiation.
Before thermal inflation, the universe is assumed to be radiation dominated. Each component evolves as where ǫ X (= ǫ X 1 + ǫ X 2 ) is the fraction of the radiation which becomes dark radiation X. For simplicity, the decoupling is assumed to occur instantaneously. Dark radiation X is divided into two components X 1 and X 2 , the former of which decays into the radiation in the visible sector. Both X 1 and X 2 contribute to the anisotropic stress, r.h.s. of (2.20).
The result of numerical calculation on the GW spectrum is shown in figure 9. Here, we consider the case where all the vacuum energy eventually goes into X particles which initially have short mean free path. Here, we used ρ r /ρ tot = 0.63 at the phase transition. We have also taken to fix t decouple and t decay , where T decay is defined as H(T = T decay ) = mΓ/T decay . Here m is the mass of X and Γ is the decay rate of X at rest. We have assumed ǫ X 1 = ǫ X 2 = ǫ X /2 and also fixed ǫ X so that ∆N eff becomes 0.5. It is seen that the spectrum changes steeply at k ≃ k PT due to the brief period of inflation. Due to the presence of dark radiation X, there is a suppression caused by the anisotropic stress for k 10 −2 k PT because GWs with such wavenumber enter the horizon after the decoupling of X. Finally, since a part of X decays into visible radiation, the suppression becomes weaker for GWs with k 10 −4 k PT which enter the horizon after the X 1 decay.

Case 2
Now, let us consider the Case 2, which has the following thermal history. (The first three events are the same as in the Case 1.) 1. A brief period of thermal inflation is caused by a scalar field φ.
2. After the phase transition, φ instantaneously decays into radiation.
3. X particles decouple from the thermal bath, after which X particles behave as dark radiation.
4. Some non-relativistic matter begins to dominate the universe.
5. The non-relativistic matter decays into the visible radiation.
Part of visible radiation or X may provide the non-relativistic matter if it becomes nonrelativistic due to the redshift. In order to study the case in which the universe evolves from the RD epoch to the MD epoch, we adopt the following equations to follow the evolution of the background: In addition, in our calculation, the massive particle which dominates the universe is assumed to originate from the visible radiation. If it is part of the non-interacting radiation we may not apply the derivation of the wave equation in ref. [28] and the calculation would be very complicated. However, from the fact that massive particles do not generate anisotropic stress, the GW spectrum in such a case is expected to be almost the same in the above figure.

Case 3
Next, let us consider the following cosmological scenario.

A brief period of thermal inflation is caused by a scalar field φ.
2. After the phase transition, φ instantaneously decays into dark radiation X.

Part of dark radiation decays into visible radiation.
Here, the decay of φ is approximated to occur instantaneously so that dark radiation has monochromatic spectrum (with energy E). Then, each component evolves as A crucial difference from the previous two cases is that φ mainly decays into dark radiation so that the effect of anisotropic stress is already significant just after the phase transition. The numerical result on the GW spectrum is shown in figure 12. We have assumed ǫ X 1 and ǫ X 2 are chosen so that ∆N eff becomes 1.3 at the decoupling and then decreases to 0.5 after the decay (ǫ X 1 = ǫ X 2 = 1/2). One finds a dip around k ≃ k PT due to the anisotropic stress caused by X, as in section 3, instead of a hill seen in the previous cases. In the low frequency limit k 10 −2 k PT , the suppression by the anisotropic stress is less efficient because X 1 does not exist when such modes enter the horizon.

Case 4
Let us consider the following cosmological scenario with late-time entropy production. The first two are the same as the Case 3.

A brief period of thermal inflation is caused by a scalar field φ.
2. After the phase transition, φ instantaneously decays into dark radiation X.
4. Non-relativistic matter decays into radiation.  Non-relativistic matter exists in many models of phase transition since the coherent oscillation of the scalar field often survives after the phase transition. In this case, we use the following set of evolution equations: ρ m + 3Hρ m = −Γρ m , (4.20) where t PT ≪ Γ −1 . The resulting GW spectra are shown in figures 13 and 14. In these figures, we have varied the parameters R rad and R coh , which are defined as R rad ≡ ρ rad ρ tot just before PT = ρ rad ρ rad + ρ vac just before PT , and with ρ coh and ρ tot being the energy density of the coherent oscillation and the total energy density, respectively. For the analysis of Case 4, we take ρ m = ρ coh . We have taken R coh = 0.1 in figure 13, and R coh = 0.01 in figure 14. The value of R rad is chosen so that, with the assumption that the vacuum energy goes only into the coherent oscillation and dark radiation, the effective neutrino number ∆N eff at the phase transition and subsequent reheating era is 2. In each case ∆N eff is diluted to 0.5 after the entropy production. Note that the features of phase transition with anisotropic stress and entropy injection appear in the figure.

Case 5
In the final example, the scalar field φ remains as a coherent oscillation after the phase transition. Thus we assume the following thermal history: 1. A brief period of thermal inflation is caused by a scalar field φ.
2. After the phase transition, φ begins a coherent oscillation, which behaves as nonrelativistic matter.
3. The coherent oscillation φ decays into dark and visible radiation.

JCAP01(2014)040
Then, the relevant evolution equations are: for t > t PT , (4.24) where B r and B X are branching ratio of φ into the radiation and dark radiation, respectively. They satisfy B r + B X = 1.
In figures 15 and 16, we show the resulting GW spectrum. In figure 15 we assumed a rather long period of thermal inflation and coherent oscillation domination. In the figure, the parameters are and B X = 0. We can see that the GW spectrum scales as k −2 if the horizon entry is during the coherent-oscillation dominated era. For the modes which experience the horizon entry twice due to the thermal inflation, the GW spectrum shows oscillatory behavior with its amplitude proportional to k −4 . (See appendix B.) We define k d as the transition frequency between these two regimes. Note that R rad parameterizes how much the initial radiation has been diluted by the short inflation, and therefore determines the ratio of the GW spectrum at k ≫ k PT to that at k = k d . On the other hand, R dil gives how much the radiation which exists at the phase transition has been diluted by the subsequent matter domination, and determines the ratio of the GW spectrum at k = k d to at k ≪ k decay ≡ aH(t = Γ −1 ). Figure 16, on the other hand, shows the GW spectrum with a short period of thermal inflation and coherent oscillation domination. Three lines (red-solid, green-dashed and bluedotted) correspond to the parameters R coh = 1, R rad = 0.67, 0.5, 0.33, R dil = 0.5, and we have taken B X so that we obtain ∆N eff = 0 (0.5) at present in the left (right) figure. Although the spectral shapes for the cases with and without dark radiation look similar, the detailed structures are different. Thus, precise observations of the height of the spectrum at k ≪ k PT and the amplitude of the oscillation around k ≃ (a few) × k PT may make it possible to tell the existence of the anisotropic stress.

Model
In this section we show that some particle-physics-motivated models can realize the Cases 1-5 discussed in the previous section. We take up two models for example, a SUSY PQ model and a SUSY majoron model.   ρ vac /ρ r =1 ρ vac /ρ r =2 Figure 16. Left: GW spectrum for Case 5 for R coh = 1, R dil = 0.5 and R rad = 0.67, 0.5, 0.33 for the red-solid, green-dashed, blue dotted lines, respectively. We also assumed ∆N eff = 0. Right: same as left but for ∆N eff = 0.5.

SUSY Peccei-Quinn model
We consider a SUSY PQ model [32]. In particular, we consider the model with the following superpotential: where Q, Q ′ andQ,Q ′ are PQ quarks which are in the fundamental and anti-fundamental representations of color SU (3), respectively, and f is the PQ scale. The superpotential has a PQ symmetry: Φ → Φe iα ,Φ →Φe −iα , Q → Qe −iα and Q ′ → Q ′ e iα . The indices i and j run -21 -JCAP01(2014)040 from 1 to a 1 and a 2 , respectively. 7 In order for the PQ symmetry to be anomalous under the color SU(3) for solving the strong CP problem, we need a 1 = a 2 . For simplicity, we assume that the Yukawa couplings y 1 and y 2 are independent of the indices i and j. The soft SUSY breaking potential is given by Here we have assumed that the soft masses m are the same for Φ,Φ and S for simplicity and neglected trilinear terms. 8 We rewrite the two superfields Φ andΦ in terms of the axion multiplet A and heavy multiplet A H as where tan θ = Φ / Φ . The axion a and saxion σ are embedded in the scalar component of A as A = (σ + ia)/ √ 2. As long as the SUSY breaking mass is much smaller than the PQ scale, Φ andΦ are constrained on ΦΦ = f 2 .
For sufficiently high cosmic temperature, the thermal effects keep Φ andΦ at the origin. As the temperature decreases, Φ andΦ begin to roll down toward the minimum. The mass matrix for Φ andΦ around the origin, including thermal effects, is given by where y ∼ √ a 1 y 1 ∼ √ a 2 y 2 . Note that, if a 1 and a 2 are larger than 1, y 1 is possible even if y 1 , y 2 1. Thus the temperature at the time of phase transition is evaluated as where, in our analysis, we assume that y λ. The ratio of radiation energy density to the total one just before the phase transition is then given by If there is no production of dark radiation, the ratio R rad also corresponds to the relative magnitude of the GW spectrum, Ω GW (k ≫ k PT )/Ω GW (k ≪ k PT ), as discussed in the previous sections.
After the phase transition, we only consider the effects of the particles in the axion multiplet A. Other particles in the PQ sector have masses larger than T PT if y 1, which we will assume in this subsection. Then, we may eliminateΦ in terms of Φ by using ΦΦ = f 2 , and obtain the potential for the scalar component of Φ. For the study of the behavior of the axion multiplet in the early universe, it is important to take account of various thermal effects [36][37][38][39]. In particular, the effective potential of Φ is given in the following form:

JCAP01(2014)040
Here, V 0 is the soft-mass term In addition, V L and V T represent thermal effects. V L is the thermal log potential [40], which comes from the fact that the masses of PQ quarks depend on the amplitude of Φ: Note that V L is effective when the temperature is lower than the masses of PQ quarks. Furthermore, V T is the thermal potential where mQ(Φ) and m Q (Φ) are masses of PQ squarks and PQ quarks, respectively. Note that V T is effective when the PQ quarks are in the thermal bath, in contrast to V L . The minimum of the potential depends on the temperature. First, in the low-temperature limit, |Φ| min ∼ f . If the thermal log potential is efficient, the minimum is given by |Φ| min ∼ α 3 T 2 /m. Finally, in the high-temperature limit, the thermal mass term determines the minimum at |Φ| min ∼ (yT f 2 /m) 1/2 . The trajectory of the PQ fields (in particular, Φ andΦ) during the phase transition depends on a i and y i . Therefore it may be the case that a considerable fraction of the initial vacuum energy may be transferred to the energy density of the saxion coherent oscillation. By using the parameter R coh given in (4.23), we phenomenologically parameterize the fraction of ρ vac transferred into the coherent oscillation energy density ρ coh ; the fraction 1 − R coh of ρ vac , which corresponds to the energy density of the oscillation of the field transverse to the trajectory ΦΦ = f 2 , is assumed to be instantly transferred into the energy densities of the axion and saxion: Hereafter we consider the case where the initial vacuum energy does not go into the visible radiation. If R coh is not negligible, the amplitudes of Φ orΦ can take large value due to the dynamics along the saxion direction. The maximal value of |Φ| during the coherent oscillation is Using m Q ′ = yf 2 / |Φ| and T PT ∼ λf /y, the PQ quarks are not in thermal equilibrium if If this condition is satisfied, the saxion may be trapped at the local minimum and the onset of the saxion coherent oscillation may be delayed.

JCAP01(2014)040
The subsequent evolution of the system depends on the decay rate and dissipation rate of the fields. First, let us consider the collision rate of the particles in the PQ sector. Without PQ quarks in the thermal bath, the axion, saxion, and axino hardly communicate with the visible sector particles. Then, they may be sequestered from the visible sector and the temperature of the PQ-sector particles may be different from that of the visible-sector particles. In particular, if the temperature of the PQ sector particles, denoted as T (PQ) , is lower than ∼ λf , the collision rate among the PQ sector particles is approximately given by Thus, when T (PQ) becomes lower than the decoupling temperature T Here the ratio of the energy density of the visible sector to the invisible one is R rad /(1 − R rad ) in the absence of the saxion coherent oscillation, therefore the temperature of the visible sector at the time of the decoupling of the PQ-sector particles is We also comment here that the effect of dissipation may be important in studying the evolution of the saxion in thermal bath [38]. When PQ quarks are not in thermal equilibrium (m Q T ), the dissipation rate of the axion, saxion and axino via the scattering with the gluon is [41] Γ diss ∼ 9α With the choices of parameters for the following discussion, however, Γ diss is smaller than the expansion rate of the universe so that the effects of dissipation are irrelevant. It means that the scattering rate between the PQ sector and visible sector is so small that the PQ sector is sequestered from the visible sector. Thus we can define the temperature of the PQ sector and the visible sector separately. After the decoupling, the saxion decays at the late epoch. The partial decay rates of the saxion into axion and gluon pairs are where µ is the higgsino mass. The decay temperature of the saxion is estimated by solving the following equation: where E σ (T ) is the averaged energy of saxion at the cosmic temperature T . The decay temperature depends on the dominant decay mode of the saxion, and also on the typical energy of the saxion at the time of the decay. If R coh ∼ 0 and R rad ∼ 0 and the saxion dominantly decays into the Higgs pair, the temperature of the visible sector T decay just after the decay is estimated as where we approximated that E σ ∼ T (PQ) and used g (PQ) * = 3.75 as the effective number of relativistic degrees of freedom in the PQ sector, and also assumed that the saxion decays while being relativistic. If the saxion decays after becoming non-relativistic, we obtain (5.25) Now let us see that the SUSY axion model can realize the Cases 1 and 2 studied in section 4.1 and 4.2, respectively. The role of φ is played by the PQ scalar, Φ andΦ. After the PQ phase transition, they decay into particles in the axion multiplet, which are thermalized due to their self interactions. (Here, we assume that the thermal bath of the PQ sector particles are sequestered from that of the visible sector.) At some point, the interactions of the PQ sector particles become very weak. Then, the mean free path of the axion and saxion becomes so long that they behave as dark radiation X. After some periods, the saxion decays into radiation. If the saxion decays while it is relativistic, it corresponds to the Case 1 and if it is non-relativistic, the Case 2 is realized.
A sample parameter set corresponding to the Case 1 is The cosmological evolution goes as follows.
• After the phase transition, the radial oscillation of Φ andΦ quickly decays into the particles in the axion multiplet.
• After the temperature drops down to T ≃ T decouple , the collision rate in the PQ sector becomes smaller than the Hubble expansion rate. So axion, saxion and axino freely stream and hence behave as dark radiation.
• If the field Φ orΦ couples to the Higgs sector, the saxion can mainly decays into the Higgs bosons. The saxion decays at T = T dec ≃ 10 6 GeV before it becomes massive. This eliminates about one-third of dark radiation. If the axino also decays without dominating the universe, ∆N eff ∼ 1 at the present universe.
The features of the GW spectrum expected with the above setup are the same as figure 9, but the positions of the characteristic wavenumber are different. In figure 9, T decouple = 10 −2 T PT and T decay = 10 −4 T PT is assumed while in the present setup they become T decouple = 10 −1 T PT and T decay = 10 −6 T PT . A sample parameter set corresponding to the Case 2 is, The cosmological evolution goes as follows.
• After the phase transition, the radial oscillation of Φ andΦ quickly decays into the PQ sector.
• After the temperature drops down to T ≃ T decouple , the collision rate in the PQ sector becomes smaller than the Hubble expansion rate, and the axion, saxion and axino become dark radiation.
• The saxion becomes non-relativistic at T ≃ m = 10 5 GeV. Then, it decays into the Higgs bosons at T dec ≃ 10 4 GeV.
The features in the GW spectrum with this setup are expected to be the same as in figures 10 and 11, although the characteristic wavenumbers are different.
If a large amount of the vacuum energy once goes into the coherent oscillation of the saxion, the Case 5 may also be realized in the SUSY PQ model. With a proper choice of -26 -

JCAP01(2014)040
parameters, most of the vacuum energy goes to the coherent oscillation of the saxion. The saxion dominates the universe and then decays; significant amount of the axion, which plays the role of dark radiation, may be produced by the decay with a relevant choice of parameters.
In this setup, the initial radiation is first diluted by the vacuum energy just like in the previous cases, and then diluted by the matter (i.e., saxion) domination which begins just after the phase transition. The former dilution can be parameterized by R rad previously defined. Let us define the quantity R dil , which parameterize the latter dilution: This determines the ratio of the spectral height at k = k d to that at k ≪ k decay ≡ aH(t = Γ −1 ). The first example in the Case 5 discussed in section 4.5 is realized with the following parameters: Thus, the cosmological evolution goes as follows: • The thermal inflation caused by the vacuum energy of Φ andΦ dilutes the initial visible radiation by a factor of 10 −4 .
• After the phase transition, almost all the vacuum energy goes to the saxion coherent oscillation. The heavy quarks remain massive since (5.14) is satisfied.
• After the coherent oscillation dilutes the radiation further by a factor of 10 −4 , it decays into visible radiation.
The second example in the Case 5 is realized with the following parameters: The cosmological evolution goes as follows: • The thermal inflation caused by the vacuum energy of Φ andΦ dilutes the initial visible radiation by a factor of O(0.1).
• After the phase transition, almost all the vacuum energy goes to the saxion coherent oscillation. The heavy quarks remain massive since (5.14) is satisfied.
• The coherent oscillation dilutes the radiation by a factor of O(0.1), and it soon decays into visible radiation.

SUSY majoron model
Next, we consider a SUSY majoron model [42,43]. The mass of the right-handed neutrino, which is needed for the seesaw mechanism [44][45][46][47], can be generated as a consequence of U(1) symmetry breaking. This U(1) may be identified as a gauged B − L symmetry, or it can be global lepton number symmetry. Here we consider the latter option. Then a Nambu-Goldstone (NG) boson appears which is associated with the spontaneous breakdown of the U(1) global symmetry. This NG boson is called majoron. (We call other fields in the majoron supermultiplet as smajoron, which is a real scalar field, and majorino, which is the fermionic superpartner of the majoron.) Let us consider the following superpotential as an example of a SUSY majoron model, 34) where N i are the right-hand neutrinos, while Φ andΦ are lepton-number violating fields. In the following argument, we take y i = y ′ i = y for simplicity. In addition, we also introduce N ′ i which have opposite lepton number. In contrast to the SUSY PQ model, this model does not lead to the thermalization of the majoron sector with the SM sector since the neutrino Yukawa coupling constants are taken to be small enough. Note also that the parameter f , the U(1) symmetry breaking scale, is not severely constrained in contrast to the PQ scale.
In the following, we assume that the majoron and the right-handed neutrinos are initially in the thermal bath. (We use the same notation as in the previous section for T PT , T decouple , T decay , R coh and R rad .) Now let us see that the SUSY majoron model can realize the Case 3 and 4 studied in The cosmological evolution goes as follows: • After the phase transition, the radial oscillation of Φ andΦ quickly decays into the majoron sector. The majoron sector does not thermalize because of the weakness of the interaction of the majoron.
• The smajoron decays at T = T dec ≃ 10 6 GeV before it becomes massive. This eliminates about one-half of the dark radiation, realizing the present ∆N eff = O(0.1).
For the Case 4, we take the following parameters: The cosmological evolution goes as follows: • After the phase transition, a small fraction of the vacuum energy goes into the coherent oscillation of the smajoron. The rest, the radial oscillation of Φ andΦ, quickly decays into the majoron sector particles. The majoron sector does not thermalize.
• The coherent oscillation of the smajoron begins to dominate the universe after the temperature decreases to 10 −1 T PT and 10 −2 T PT in figures 13 and 14, respectively. This dilutes the dark radiation produced at the phase transition. Note that the smajoron particles are still relativistic or just becoming non-relativistic at these temperatures.
• The coherent oscillation (and smajoron particles in the case of figure 14), decays into visible radiation and the present ∆N eff becomes O(0.1).

Conclusions
In this paper we have studied various effects of the cosmological events in the early universe on the inflationary GW spectrum. In particular, some (well-motivated) models of cosmological phase transition, including the PQ phase transition, in general lead to complicated structures in the GW spectrum, which might be detectable in future space-based GW detectors [15][16][17][18].
Hence, if the spectrum of the inflationary GWs is determined by future experiments, it will give clues to the high-energy physics beyond the reach of accelerator experiments.

32πG
h ij;0 h ij;0 osc (t, x). (A.11) We define the effective energy of GWs per logarithmic frequency ρ GW (t, k) through ρ GW (t) ≡ d ln k ρ GW (t, k), (A.12) where we have assumed the homogeneity and isotropy of the primordial GWs. Because of the subhorizon condition k ≫ aH, we may neglect terms with time derivative on the scale factor in (A.11). If the universe is filled with perfect fluid at t, (2.7) tells us that h is a harmonic oscillator with frequency k/a. Thus for subhorizon modes, we can use the relation ḣ 2 ij osc = (k/a) 2 h 2 ij osc to write ρ GW (t, k) as ρ GW (t, k) = 1 32πG k 2 a 2 k 3 2π 2 P h (t, k) = 1 32πG k 2 a 2 ∆ 2 GW (t, k). (A.13) 10 Here we have considered GWs entering the horizon at the RD era. For those entering in the MD era, the factor 1/2 in eq. (A.7) should be replaced with 9/32, which is derived from the analytical solution given in eq. (2.11).

JCAP01(2014)040
The GW spectrum Ω GW (t, k) is defined as the fraction of ρ GW (t, k) to the critical density: Ω GW (t, k) ≡ ρ GW (t, k) ρ tot (t) . (A.14) The present value of Ω GW may be decomposed as follows Assuming that no significant entropy injection has occurred between the time of horizon entry for GWs with wavenumber k and the matter-radiation equality, we obtain