[type=editor, auid=000,bioid=1] [1] [2] [1]

[cor1]Corresponding author

# Fast FPGA algorithm for neutron-gamma discrimination

Haoqi Ye  Lian Chen  Ge Jin  State Key Laboratory of Particle Detection and
Electronics, University of Science and Technology of China, Hefei 230026,
China<EMAIL_ADDRESS>
###### Abstract

Various pulse shape discrimination methods have been used to solve the
neutron-gamma discrimination problem. But most of them are limited to off-line
calculation due to the computation amount and FPGA performance. In order to
realize real time discriminating neutron and gamma, a new algorithm based on
the traditional pulse shape discrimination methods was proposed in this paper.
The new algorithm takes into account the physical properties of the pulse
signal, which greatly reduces the computation and dead time without losing the
precision, and can work on FPGA directly. It has a good performance in the
actual experiment based on CLLB scintillation detector.

###### keywords:

neutron-gamma discriminationonline PSD method FPGA

##  1 Introduction

Neutron sources are always accompanied by the de-excitation gamma rays whose
discrimination from neutrons is basically a complicated task[1]. To solve this
problem, in the liquid scintillator and plastic scintillator detectors, which
have been commonly used in neutron detection recently, the pulse shape
discrimination(PSD) technique is used for n/$\gamma$ ray resolution. In such
applications, PSD technique is used to separate neutron interaction events
from interfering $\gamma$-ray interactions by exploiting a difference in the
intensity of the slow component of the light pulses in organic scintillators
that are initiated by the recoil protons and electrons. In the past few
decades, various PSD techniques have been developed. The earliest was PSD
circuit. With the development of digital circuits, in recent years, there have
been a variety of new PSD methods based on digital circuits and computer
technology. such as pattern reorganization method (PRM)[2], wavelet packet
transform method (WPTM) [3] and discrete Fourier transform method(DFTM) [4].
However, due to the resource limitations of ASIC and FPGA, most of the PSD
method are implemented in computers. In most of the traditional pulse-shape
discrimination, the data is collected by data acquisition cards, and after
being preprocessed by FPGA, packaged and uploaded to the computer, then the
data is calculated offline. These algorithms have problems such as wasting
FPGA performance and increasing data processing time. Therefore the new PSD
method is needed which is suitable for FPGA characteristics. This PSD method
needs small computation, high speed, high performance, and should work online
on FPGA .

##  2 Traditional PSD method

###  2.1 Charge-comparison method

In scintillators, the optical pulses excited by charged particles have
different attenuation time components. The intensity ratio of the fast and
slow components is related to the mass and charge of the excited particles,
that is, to the ionization density formed by charged particles in the
scintillator. For example, in organic scintillator stilbene crystal, the
intensity ratio of the slow component to the fast component in the light pulse
excited by fast neutrons is 4.5 times that of the light pulse excited by
$\gamma$ rays[5]. The fluorescence emitted by the scintillator is collected by
a photomultiplier tube. When the photomultiplier tube operates in a linear
range, the shape of the current pulse drawn from its anode reflects the shape
of the light pulse emitted by the scintillator, i.e., the current pulse can be
expressed as (1):

$I\left(t\right)=I_{f}\left(\rho\right)e^{-\frac{t}{\tau_{f}}}+I_{s}\left(\rho\right)e^{-\frac{t}{\tau_{s}}}$
(1)

In equation, $\tau_{f}$,$\tau_{s}$,$I_{f}(\rho)$,$I_{s}(\rho)$ are the decay
time of fast and slow components and the maximum value of current pulses. The
total charge can be obtained by integrating the current pulse with an integral
loop with a large time constant

Q=∫_0^∞I( t ) dt=I_f( ρ)τ_f+I_s( ρ) τ_s

=Q_f( ρ) +Q_s( ρ)

(2)

The charge pulse integration is also composed of two parts: fast and slow.
Neutrons and gamma rays form different $Q_{s}/Q_{total}$, so this ratio is
used to identify particles. This method is called charge comparison method. In
actual calculation, the parameter PSD is defined as the ratio of the integral
of the trailing edge of the pulse to the integral of the whole pulse. To
quantify the n/$\gamma$ discrimination power at a given energy threshold, a
figure of merit(FOM)is used. The _FOM_ parameter is defined as (3):

$FOM=\frac{\tau}{FWHM_{n}+FWHM_{\gamma}}$ (3)

The _FOM_ represents the performance of event discrimination rate. The
shortcoming of this algorithm is that it requires a lot of computation and
takes up a lot of resources, which is not conducive to online identification
on FPGA platform.

Figure 1: The PSD performance of the system with the charge-comparison method.
Figure 2:  Schematic diagram of PSD spectrum.FWHMγ and FWHMn are half-widths
of two PSD peaks.$\tau$ is the distance between two peaks. Figure 3: The
voltage signal obtained after processing the current signal of the neutron and
gamma rays at the output of the photomultiplier tube. The charge comparison
method uses the ratio of the integral from $t_{1}$ to $t_{2}$ of the signal
waveform to the integral from $t_{0}$ to $t_{2}$ to distinguish neutrons and
gamma rays

###  2.2 Pulse gradient analysis method

Pulse gradient analysis (PGA) method is a simple discriminating method based
on the different shapes of neutron and $\gamma$-ray signals[6]. This method
selects the pulse peak and a sample point after the peak to calculate the
pulse gradient. In general, the optimal value range of the time span T between
the pulse wave peak and the selected sample point is 15-25ns.The specific
value depends on the material of the scintillation detector and the
characteristics of the photomultiplier tube. This method has low computation
and high speed, but it is very sensitive to the random fluctuation of the
signal due to the great influence of noise

##  3 Partial Charge-to-Peak ratio method

In the charge-comparison method, the ratio of the back integral to the total
integral of the signal was used in the identification of the fast and slow
components. When calculating the PSD value, the denominator in the formula (1)
plays a role of normalization. However, as the denominator is the full
integral of voltage, which requires more resources, we need a substitute
physical quantity. Considering the physical meaning of the pulse waveform, the
full integral at the denominator is the integral of the voltage signal
generated from the operational amplifier by the pulse output from the
photomultiplier tube. We can consider it as a value that is proportional to
the magnitude of the charge (that is, the integral of the current waveform).
From the reaction principle of neutrons and gamma rays with scintillator [5],
the waveform of the normalized output current of photomultiplier tube is only
related to the type of reaction. Therefore, for a definite reaction, the
voltage pulse integral is the value proportional to the peak of the current
waveform In the actual calculation, $V_{peak}$, which is algebraically related
to $I_{peak}$, is used. The relation between $V_{peak}$ and charge integral is
(4)

$\int_{t_{0}}^{t_{2}}{Vdt=\alpha V_{peak}}$ (4)

The definitions of $t_{0}$,$t_{2}$ are shown in Fig.2. In this condition, the
peak voltage is used instead of the charge integral. This not only reduces the
amount of calculation, but also accords with the physics principle during
calculation. According to the theory of fast and slow components, the fast
component has a greater influence near the peak value, and the signal-to-noise
ratio is poor at the tail of the pulse. So for the back integral, we don’t
have to take into account the integration of both ends when we integrate the
pulse, which means that the equation becomes the integral from
$t_{1}^{\prime}$ to $t_{2}^{\prime}$ as shown in Fig 4

Figure 4: The voltage signal obtained after processing the current signal of
the neutron and gamma rays at the output of the photomultiplier tube

The advantage of this modification is that for the pulse with a long trailing
edge, the calculation can be greatly reduced and a lot of storage resources
can be saved by selecting $t_{1}^{\prime}$ and $t_{2}^{\prime}$. The ratio of
the back integral of the neutron/gamma waveforms before and after modification
are (5) and (6):

$\eta=\frac{\int_{t_{1}}^{t_{2}}{V_{N}dt}}{\int_{t_{1}}^{t_{2}}{V_{\gamma}dt}}$
(5)
$\eta^{\prime}=\frac{\int_{t_{1}^{\prime}}^{t_{2}^{\prime}}{V_{N}dt}}{\int_{t_{1}^{\prime}}^{t_{2}^{\prime}}{V_{\gamma}dt}}$
(6)

Since the influence of some fast components is eliminated, the modified back-
edge integral ratio is lower than the original back-edge integral
ratio($\eta^{\prime}>\eta$). The influence brought by noise is reduced
considerably and the signal integral is basically unchanged. Based on the
physical principle of the reaction, the mathematical relationship between the
physical quantities and the requirement of the calculation formula for FPGA,
the PSD calculation formula of Partial charge-to-Peak ratio method on the
basis of charge-comparison method and pulse gradient analysis method is shown
as (7):

$PSD=\frac{\int_{t_{1}^{\prime}}^{t_{2}^{\prime}}{Vdt}}{V_{peak}}$ (7)

Compared with the algorithm of the original PSD method, by setting the values
of $t_{1}^{\prime}$ and $t_{2}^{\prime}$, 80$\%$-90$\%$ of the computation and
70$\%$ of the required storage resources can be reduced. In order to study
discrimination performance, we consider a pair of gamma rays and neutrons with
the same peak current pulse. The ratios of the neutron and gamma pulse PSD
values of the old and new methods are (8) and (9):

$\frac{\left(PSD_{new}\right)_{N}}{\left(PSD_{new}\right)_{\gamma}}=\frac{\frac{\int_{t_{1}^{\prime}}^{t_{2}^{\prime}}{V_{N}dt}}{V_{peak}}}{\frac{\int_{t_{1}^{\prime}}^{t_{2}^{\prime}}{V_{\gamma}dt}}{V_{peak}}}=\frac{\int_{t_{1}^{\prime}}^{t_{2}^{\prime}}{V_{N}dt}}{\int_{t_{1}^{\prime}}^{t_{2}^{\prime}}{V\gamma
dt}}=\eta^{\prime}$ (8)
$\frac{\left(PSD_{old}\right)_{N}}{\left(PSD_{old}\right)\gamma}=\frac{\frac{\int_{t_{1}}^{t_{2}}{V_{N}dt}}{\int_{t_{0}}^{t_{2}}{V_{N}dt}}}{\frac{\int_{t_{1}}^{t_{2}}{V\gamma
dt}}{\int_{t_{0}}^{t_{2}}{V\gamma
dt}}}=\frac{\frac{\int_{t_{1}}^{t_{2}}{V_{N}dt}}{\alpha_{N}V_{peak}}}{\frac{\int_{t_{1}}^{t_{2}}{V\gamma
dt}}{\alpha_{\gamma}V_{peak}}}=\frac{\alpha_{\gamma}}{\alpha_{N}}\eta$ (9)

From Fig.1, the area under $\gamma$ pulse waveform is smaller than that under
neutron pulse when $V_{peak}$ is same, so $\alpha_{N}$>$\alpha_{\gamma}$ , and
$\eta^{\prime}>\eta$ has been known previously, and thus we can obtain (10).

$\frac{\left(PSD_{new}\right)_{N}}{\left(PSD_{new}\right)_{\gamma}}\ >\
\frac{\left(PSD_{old}\right)_{N}}{\left(PSD_{old}\right)_{\gamma}}$ (10)

In addition, this algorithm reduces the influence of noise. Theoretically, the
PSD performance of Partial charge-to-peak ratio method should be better than
that of charge-comparison method. Compared with the pulse gradient method, the
calculation formulas of the two methods are similar. However, considering the
noise, random error and other factors, the Partial charge-to-peak ratio method
is more stable than the pulse gradient method under the condition that the
calculation amount is not too different.

##  4 FPGA Programs

###  4.1 The optimization of Algorithm

In FPGA, we express the waveform integral as the sum of the values of the
discrete points from ADC. Because the integral position depends on the peak
position, we need to find the peak before we integrate the pulse. In this
case, the pulse signal should be preserved first which requires a lot of FPGA
resources and a long dead time for integration. Three registers used to store
the peak, the difference between current time and peak time, and the integral
value are used to solve this problem. In order to reduce the dead time of the
system, integration is carried out at the same time as receiving pulse
signals. Setting the value of ADC input to FPGA as $V_{i}$, when $V_{i}$
exceeds the set threshold $V_{T}$, the program starts to work. We use a
register to record the time difference between the current signal value and
the peak value. When the time difference is between $T_{1}^{\prime}$ and
$T_{2}^{\prime}$, the current signal value is added to the value in the
register used to record the integral. If the current signal value exceeds the
peak value, the values in register for time difference and register for
integral are cleared to zero, and the peak is updated to the current value.
When the value in the time register exceeds the preset pulse length $T_{e}$,
the PSD parameter value should be calculated. The PSD parameter value is
defined as integral divided by peak. The $V_{T}$, $T_{1}^{\prime}$,
$T_{2}^{\prime}$, $T_{e}$ are defined as Fig.5 shows.

Figure 5: Discrete pulse signal which FPGA gets from ADC. $T_{1}^{\prime}$ and
$T_{2}^{\prime}$ are the upper and lower limits of the integration interval,
and $T_{e}$ is the pulse length, usually set as the time at which the signal
value is reduced to 5$\%$ of the peak value. i is the time difference between
the current time and the peak time, $V_{peak}$ is the current peak value.

Under this optimization, it is not necessary to record the whole waveform
signal. Only a few register are needed to record the time difference, signal
maximum value and integral. Considering the requirements of online work, only
two judgment statements are considered per clock time, and only one or two
calculations are performed to minimize the dead time. This design can not only
save a lot of resources, but also perfectly accord with the characteristics of
small amount of computation and short dead time required by FPGA.

###  4.2 FPGA state machine

According to the algorithm, the FPGA state machine is shown in Fig6

Figure 6: FPGA state machine, $S_{0}$-$S_{5}$ are state codes, $V_{i}$ is the
current signal value, and Q is the integral. Other parameters are the same as
in Fig.5.

##  5 Experimental results and discussion

We used the experimental platform as Fig.7 shows. This design uses a
1.5inch*1.5inch CLLB crystal packaged with a R9420 PMT made by Hamamatsu
Company. The readout board reads the charge signal from photocathode, and then
transfer the signal to Front End Board (FEB). After being shaped and amplified
by FEB, the signal is digitized by Data Collect and Process (DCP) board.

Figure 7: Experimental platform based on CLLB crystal packaged with PMT

A 250MSPS,14bit ADC in DCP board is used to digitize the signal pulse, then
online PSD methods deployed by a FPGA process the data stream. At the same
time, FPGA upload the data stream to computer. We use other offline PSD method
to process the data in computer. In the FEB, in order to shape the original
waveform ,we used a SK filter to spin wave to suit our 250MSPS ADC. In this
case, the trailing edge of the pulse is about 1000ns(95$\%$-5$\%$)because the
ADC sampling rate is not high. To analyze the properties of different PSD
methods, we used Partial Charge-to-peak ratio method and pulse gradient
analysis method as online methods, used charge-comparison method as offline
method. For the $\varDelta t$ parameter of gradient analysis method, the most
appropriate value 600ns, which is 150 clock in ADC, was selected after many
experiments. For the $T_{1}^{\prime},T_{2}^{\prime}$ of partial charge-to-peak
ratio method, we chose $T_{1}^{\prime}$=150(600ns),
$T_{2}^{\prime}$=180(720ns), $\varDelta T$=30(120ns). The energy threshold was
chosen the same value (300keV). Considering the balance of precision and
resources, we selected 1024 channels precision in FPGA, that is, there were
1024 channels for storing PSD parameter. The value of the divided charge
integral was shifted to the left before the division to accommodate the
multichannel. The PSD spectrums of these three PSD methods are shown in
Fig.8-10,The _FOM_ parameter was determined by fitting a double Gaussian
function to the distribution of the PSD parameter.

Due to the lead shielding between the detector and the neutron source in the
experiment, the detector received a very small dose of gamma rays, resulting
in the _FOM_ of the PSD energy spectrum in Gaussian fitting much smaller than
that without lead shielding As can be seen from the results in Fig.6-8,
partial charge-to-peak ratio method has better performance compared with
gradient analysis method and charge-comparison method. In addition, it has a
small amount of computation and adapts to the characteristics of FPGA
platform, which indicates that this method is of great practical significance.
Another study of this method is to consider the values of parameters. We need
to consider energy threshold, the place where start the integration
$T_{1}^{\prime}$ and the place where end the integration $T_{2}^{\prime}$.
From the noise point of view,the increase in the _FOM_ with increasing energy
threshold is apparent. About the influence of other parameters, the larger
$T_{1}^{\prime}$ is, the smaller the effect of the fast component is; the
smaller $T_{2}^{\prime}$ is, the smaller the effect of noise is, larger the
$T_{2}^{\prime}-T_{1}^{\prime}$ value is, more stable the performance is (the
smaller the influence of random error) but also the greater the computational
effort is. In order to study the influence of each parameter on the _FOM_ , we
made a series of experiments and calculations. Results are shown in Table 1.

Figure 8: Charge-comparison method’s PSD spectrum with _FOM_ , the PSD size has been adjusted. Figure 9: Gradient analysis method’s PSD spectrum with _FOM_ Figure 10: Partial charge-to-current ratio method’s PSD spectrum with _FOM_ Table 1: _FOM_ s obtained by selecting different $T_{1}^{\prime},T_{2}^{\prime}$ with threshold =300keV .The $\varDelta T$ in Table1 means $T_{2}^{\prime}-T_{1}^{\prime}$ .The number means how many clocks. 1 clock = 4ns ($F_{clk}$ = 250Mhz) | $\Delta T=10$ | $\Delta T=30$ | $\Delta T=50$
---|---|---|---
$T_{1}^{\prime}=100$ | 0.75 | 0.84 | 0.91
$T_{1}^{\prime}=150$ | 1.02 | 1.06 | 1.03
$T_{1}^{\prime}=200$ | 0.9 | 1.02 | 1.00

The selection of appropriate integral interval has a great influence on the
_FOM_. In this experimental platform, the appropriate integral interval is
600ns-720ns. In other experimental devices, the appropriate interval should be
selected according to the detector and the forming function of the forming
circuit. It is appropriate to select 10$\%$ pulse length at about 60$\%$ of
the pulse falling edge (95$\%$-5$\%$). The experimental results of the energy
threshold value are shown in Table 2.

Table 2: _FOM_ s under different thresholds with $T_{1}^{\prime}$ = 150(600ns),$\varDelta T$ = 30(120ns) Threshold | 300keV | 400keV | 480keV
---|---|---|---
$T_{1}^{\prime}=150,\varDelta T=30$ | 1.06 | 1.12 | 1.16

As we expected, the _FOM_ increased as the threshold increased. The selection
of threshold value should be made according to the specific situation.
Excessive threshold value will bring more data loss. One limitation of FPGA is
that it cannot use floating numbers. So the same algorithm has different
accuracy on FPGAs and computers. When calculating the PSD value, the integral
is shifted to the left first. This processing is equivalent to the precision
setting of floating numbers. The algorithm was experimented on both FPGA and
computer to compare the precise difference. The results in Table 3 show that
there is little difference between the results obtained by computer and those
obtained by FPGA.

Table 3: _FOM_ s with $\varDelta T$=30(120ns), threshold = 300keV | Computer | FPGA
---|---|---
$T_{1}^{\prime}=100$ | 0.86 | 0.84 |
$T_{1}^{\prime}=150$ | 1.07 | 1.06 |
$T_{1}^{\prime}=200$ | 1.02 | 1.02 |

##  6 Conclusions

In this study, we discuss a fast n/$\gamma$ discrimination algorithm which can
work on FPGA in real time with uncomplicated computation. By considering the
characteristics of the fast and slow components in the pulse signal and
combining the advantages of the charge comparison method and the pulse
gradient method, the calculation amount were reduced 80$\%$ by this method
compared with the former method. The method was tested with a CLLB crystal
packaged with a R9420 PMT, and generated a good PSD spectrum online in FPGA.
In the PSD resolution performance, partial charge-to-peak ratio method is
better than traditional algorithms. In parameter setting, 60$\%$-70$\%$ of the
falling edge is recommended as an integral interval. The algorithm gets
similar results running in FPGA and computer, but has little dead time(three
clocks) in FPGA, which means this algorithm is suitable for online n/$\gamma$
discrimination.

## Acknowledgments

This work was supported by the National Natural Science Foundation of China
under Grant(No. 11875249).

## References

  * [1] G. F. Knoll, _Radiation detection and measurement_. John Wiley & Sons, 2010.
  * [2] D. Takaku, T. Oishi, M. Baba _et al._ , “Development of neutron-gamma discrimination technique using pattern-recognition method with digital signal processing,” _Progress in Nuclear Science and Technology_ , vol. 1, pp. 210–213, 2011.
  * [3] M. Safari, F. A. Davani, H. Afarideh, S. Jamili, and E. Bayat, “Discrete fourier transform method for discrimination of digital scintillation pulses in mixed neutron-gamma fields,” _IEEE Transactions on Nuclear Science_ , vol. 63, no. 1, pp. 325–332, 2016.
  * [4] M. Nakhostin, “Digital discrimination of neutrons and $\gamma$-rays in liquid scintillation detectors by using low sampling frequency adcs,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 916, pp. 66–70, 2019\. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0168900218315675
  * [5] C. Yu, “Charge comparison method used to discriminate photons and neutrons in liquid scintillators(chinese with english abstract),” _High energy physics and nuclear physics_ , vol. 23, no. 7, pp. 616–621, 1999.
  * [6] B. D’Mellow, M. Aspinall, R. O. Mackin, M. J. Joyce, and A. J. Peyton, “Digital discrimination of neutrons and $\gamma$-rays in liquid scintillators using pulse gradient analysis,” _Nuclear Inst & Methods in Physics Research A_, vol. 578, no. 1, pp. 191–197, 2007.


# Modeling Cosmological Perturbations of Thermal Inflation

Jeong-Myeong Bae  Center for Theoretical Physics, Department of Physics and
Astronomy, Seoul National University, Seoul 08826, Republic of Korea  School
of Undergraduate Studies, College of Transdisciplinary Studies, DGIST, Daegu
42988, Republic of Korea  Sungwook E. Hong  Korea Astronomy and Space Science
Institute, Daedeok-daero 776, Yuseong-gu, Daejeon 34055, Republic of Korea
Astronomy Campus, University of Science and Technology, Daedeok-daero 776,
Yuseong-gu, Daejeon 34055, Republic of Korea  Heeseung Zoe
<EMAIL_ADDRESS>Department of Physics, Izmir Institute of Technology,
Gulbace, Urla 35430, Izmir, Turkiye  Department of Physics Education, Pusan
National University, Busan 46241, Republic of Korea  Department of Physics
Engineering, Istanbul Technical University, Maslak 34469, Istanbul, Turkiye

###### Abstract

We consider a simple system consisting of matter, radiation and vacuum
components to model the impact of thermal inflation on the evolution of
primordial perturbations. The vacuum energy magnifies the primordial modes
entering the horizon before its domination, making them potentially
observable, and the resulting transfer function reflects the phase changes and
energy contents. To determine the transfer function, we follow the curvature
perturbation from well outside the horizon during radiation domination to well
outside the horizon during vacuum domination and evaluate it on a constant
radiation density hypersurface, as is appropriate for the case of thermal
inflation. The shape of the transfer function is determined by the ratio of
vacuum energy to radiation at matter-radiation equality, which we denote by
$\upsilon$, and has two characteristic scales, $k_{\rm a}$ and $k_{\rm b}$,
corresponding to the horizon sizes at matter radiation equality and the
beginning of the inflation, respectively. If $\upsilon\ll 1$, the universe
experiences radiation, matter and vacuum domination eras and the transfer
function is flat for $k\ll k_{\rm b}$, oscillates with amplitude $1/5$ for
$k_{\rm b}\ll k\ll k_{\rm a}$ and oscillates with amplitude $1$ for $k\gg
k_{\rm a}$. For $\upsilon\gg 1$, the matter domination era disappears, and the
transfer function reduces to being flat for $k\ll k_{\rm b}$ and oscillating
with amplitude $1$ for $k\gg k_{\rm b}$.

##  I Introduction

Inflation provides a theoretical ground for our understanding of the universe
[1, 2, 3, 4, 5]. It makes the universe homogeneous, isotropic, and flat and
dilutes unwanted or unobserved relics such as monopoles. It has been
constrained by observations of the large-scale structure (LSS) and cosmic
microwave background (CMB) [6]. In supersymmetric cosmology, however, the
moduli fields are dangerous to the big bang nucleosynthesis (BBN), if not
effectively removed [7, 8, 9]. There have been attempts to solve the moduli
problem by arranging low-energy inflation after the primordial inflation [10],
but it is tricky to control the moduli density in the permissible range
because the moduli are regenerated after the low-energy inflation.

This moduli problem can be solved by introducing a thermal inflation [11, 12],
a brief and secondary inflationary phase after the primordial inflation, being
realized by thermal effects on flat directions in supersymmetric theories [13,
14, 15, 16, 17, 18]. The thermal inflation occurs at the primordial dark
period between the end of primordial inflation and BBN. Some mechanisms,
including reheating/preheating and baryogenesis/leptogenesis, have been
suggested for this period. Thermal inflation scenario has a very different
post-inflationary history from the standard scenario and provides new
predictions about the primordial dark period.

While there are some prospects from gravitational waves [19, 20, 21] and
collider physics [22], it is hard to directly probe these post-inflationary
mechanisms in “both” standard and thermal inflation scenarios by observations
up to now. Thermal inflation, however, gives us a better chance: it magnifies
modes that have entered the horizon during the primordial dark period and
becomes sensitive to the physics there. In [23], we studied these effects at
scales smaller than the horizon size at the beginning of thermal inflation.
Thermal inflation suppresses the power spectrum of those modes and hence gives
the suppression of CMB $\mu$-distortions [24], the 21-cm hydrogen power
spectrum at or before the epoch of reionization, and the formation of galaxy
substructures [25].

In this paper, we consider a system consisting of matter, radiation and vacuum
components to model the impact of thermal inflation on the evolution of
primordial perturbations. By assuming that the primordial inflation generates
cosmological perturbations in a standard way, we study the growth of the
perturbations affected by the secondary thermal inflation. We demonstrate how
the cosmological scenario between the primordial and thermal inflations
affects the power spectrum and possibly leaves observable small-scale
features. We consider the ratio of vacuum energy to radiation at matter-
radiation equality as a key parameter for characterizing the density
perturbation, which explains the previous result of [23]. The detail of moduli
dynamics during its domination between the two inflations also impacts the
power spectrum, but it is beyond the scope of the present paper.

This paper is organized as follows. In section II, we review the thermal
inflation scenario and introduce the setting for the perturbation analysis. In
section III, we examine a simple system of vacuum energy and radiation to
model thermal inflation. In section IV, we study a system of vacuum energy,
radiation, and matter to calculate the curvature perturbations for thermal
inflation analytically. In section V, we present numerical calculations of the
power spectrum according to the initial ratio between matter and radiation. In
section VI, we summarize the results and discuss the future work.

##  II Review of Thermal Inflation

Figure 1: The cosmological history of thermal inflation and four
characteristic scales.

In this section, we briefly review the thermal inflation scenario [11, 12] and
introduce density perturbations as a setting for the next sections (for the
details, [23]). FIG. 1 shows the cosmological history with thermal inflation,
and the universe experiences the following four major phases.

  * •

Primordial inflation: At $t<t_{\rm a}$, the primordial inflation generates the
scale-invariant power spectrum, and unknown post-inflationary era including
the reheating process follows.

  * •

Moduli domination: In supersymmetric theories, the moduli mass in vacuum is an
order of the soft supersymmetry breaking scale, i.e., $m_{\mathrm{mod}}\sim
m_{\mathrm{s}}\sim 10^{3}~{}\mathrm{to}~{}10^{4}\mathinner{\mathrm{GeV}}$. At
$t\sim t_{\mathrm{a}}$, the Hubble scale gets $H\sim m_{\mathrm{s}}$, and the
moduli start dominating over the universe. At $t>t_{\mathrm{a}}$ or
$H<m_{\mathrm{s}}$, the moduli start oscillating with Planckian amplitude, and
the oscillation may be long-lived to spoil the BBN.

  * •

Thermal inflation: At $t=t_{\rm b}$, thermal inflation begins to resolve the
moduli problem. Thermal inflation is realized by flaton(s) that is trapped by
finite temperature potential induced by radiation in the universe. The flaton
potential is
$\mathinner{V\mathopen{\left(\phi,T\right)}}=V_{0}+\frac{1}{2}\left(\sigma^{2}T^{2}-m_{\phi}^{2}\phi^{2}\right)+\cdots$
where the thermal coupling $\sigma$ is not small and the flaton mass is
$m_{\phi}\sim m_{\mathrm{s}}$. For $V^{1/4}_{0}\sim
10^{5}~{}\mathrm{to}~{}10^{8}\mathinner{\mathrm{GeV}}$, thermal inflation has
e-folds $N_{\rm bc}\sim 10~{}\mathrm{to}~{}15$ and dilute the pre-existing
moduli by $\Delta_{\rm TI}\sim e^{3N_{\rm bc}}$.

  * •

Flaton domination: As the temperature of the universe drops, thermal inflation
ends until $t\sim t_{\mathrm{c}}$, and the flaton starts oscillation around
its vacuum expectation value $\phi_{\rm vev}\sim\sqrt{V_{0}}/m_{\phi}$ and
dominates over the universe as a matter phase at $t\gtrsim t_{\mathrm{c}}$.
The moduli can be regenerated after thermal inflation but are diluted further
by $\Delta_{\rm flaton}\sim e^{N_{\rm cd}}$ due to the flaton decay.

  * •

Radiation domination: At $t\sim t_{\mathrm{d}}$, flaton decays to yield the
standard radiation dominated universe, at the temperature $T_{\mathrm{d}}\sim
10^{-2}~{}\mathrm{to}~{}10^{2}\mathinner{\mathrm{GeV}}$. The universe steps
into the standard cosmic history after these radiation domination are
recovered.

In summary, three distinct extra eras of thermal inflation - moduli
domination, thermal inflation, and flaton domination - are inserted between
the primordial inflation and radiation domination of the standard scenario.
Therefore, the cosmological perturbation evolves differently from the standard
inflation scenario.

We introduce four characteristic scales, $k_{\mathrm{a}}$, $k_{\mathrm{b}}$,
$k_{\mathrm{c}}$ and $k_{\mathrm{d}}$, where

$k_{x}\equiv a_{x}H_{x}$ (1)

corresponds to the comoving scale of the horizon at the era boundary $t_{x}$.
In FIG. 1, modes with $k<k_{\mathrm{b}}$ remain outside the horizon throughout
the thermal inflation eras and are not affected by thermal inflation. Modes
with $k_{\mathrm{b}}<k<k_{\mathrm{a}}$ enter the horizon during moduli
domination, and their growth is modified. Modes with $k>k_{\mathrm{a}}$ enter
the horizon before moduli domination and so probe that unknown era. Modes with
$k>k_{\mathrm{d}}$ reenter the horizon during flaton domination and so will be
twice modified. Modes with $k>k_{\rm c}$ never exit the horizon throughout
thermal inflation eras.

Now we review the perturbation analysis in [23]. For $t\lesssim t_{\rm a}$,
the post-inflationary physics of the reheating and modulogenesis before the
moduli domination is unknown. However, for

$\displaystyle t_{\mathrm{a}}\ll t<t_{\mathrm{c}}~{},$ (2) $\displaystyle k\ll
k_{\mathrm{a}}~{},$ (3)

we can study the perturbations based on a simple system of moduli matter (m),
thermal radiation (r) and vacuum energy ($V_{0}$). The energy density and
pressure are

$\displaystyle\rho=\rho_{\mathrm{m}}+\rho_{\mathrm{r}}+V_{0}$ (4)
$\displaystyle p=\frac{1}{3}\rho_{\mathrm{r}}-V_{0}~{}.$ (5)

The scalar part of the metric perturbation is [26]

$\tilde{ds}^{2}=(1+2A)\mathinner{dt}^{2}-2B_{,i}\mathinner{dt}\mathinner{dx^{i}}-\left[(1+2\mathcal{R})\mathinner{a^{2}\mathopen{\left(t\right)}}\delta_{ij}+2C_{,ij}\right]\mathinner{dx^{i}}\mathinner{dx^{j}}~{}.$
(6)

We describe the perturbations in moduli and radiation by introducing the gauge
invariant variables

$\displaystyle\mathcal{R}_{\delta\rho_{\mathrm{m}}}\equiv\mathcal{R}-\frac{H}{\dot{\rho}_{\mathrm{m}}}\delta\rho_{\mathrm{m}}$
(7)
$\displaystyle\mathcal{R}_{\delta\rho_{\mathrm{r}}}\equiv\mathcal{R}-\frac{H}{\dot{\rho}_{\mathrm{r}}}\delta\rho_{\mathrm{r}}$
(8)

where $\mathcal{R}_{\delta\rho_{\mathrm{m}}}$ is the curvature perturbation on
constant moduli density hypersurfaces and
$\mathcal{R}_{\delta\rho_{\mathrm{r}}}$ is the curvature perturbation on
constant radiation density hypersurfaces. We get two coupled equations for
scalar perturbations [23]

$\displaystyle\ddot{\mathcal{R}}_{\delta\rho_{\mathrm{m}}}+H\left(2+\frac{\rho_{\mathrm{m}}}{\rho+p+\frac{2}{3}q^{2}}\right)\dot{\mathcal{R}}_{\delta\rho_{\mathrm{m}}}-\frac{1}{3}q^{2}\left(\frac{\rho_{\mathrm{m}}}{\rho+p+\frac{2}{3}q^{2}}\right)\mathcal{R}_{\delta\rho_{\mathrm{m}}}$
$\displaystyle{}=-\frac{\frac{4}{3}\rho_{\mathrm{r}}}{\rho+p+\frac{2}{3}q^{2}}\left(H\dot{\mathcal{R}}_{\delta\rho_{\mathrm{r}}}-\frac{1}{3}q^{2}\mathcal{R}_{\delta\rho_{\mathrm{r}}}\right),$
(9)
$\displaystyle\ddot{\mathcal{R}}_{\delta\rho_{\mathrm{r}}}+H\left(1+\frac{\frac{8}{3}\rho_{\mathrm{r}}}{\rho+p+\frac{2}{3}q^{2}}\right)\dot{\mathcal{R}}_{\delta\rho_{\mathrm{r}}}+\frac{1}{3}q^{2}\left(1-\frac{\frac{8}{3}\rho_{\mathrm{r}}}{\rho+p+\frac{2}{3}q^{2}}\right)\mathcal{R}_{\delta\rho_{\mathrm{r}}}$
$\displaystyle{}=-\frac{2\rho_{\mathrm{m}}}{\rho+p+\frac{2}{3}q^{2}}\left(H\dot{\mathcal{R}}_{\delta\rho_{\mathrm{m}}}-\frac{1}{3}q^{2}\mathcal{R}_{\delta\rho_{\mathrm{m}}}\right),$
(10)

where $q\equiv k/a$. For $t_{\mathrm{a}}\ll t$, we have

$\rho_{\mathrm{r}}\ll\rho_{\mathrm{m}}$ (11)

and simplify Eqs. (II) and (II) to

$\displaystyle\ddot{\mathcal{R}}_{\delta\rho_{\mathrm{m}}}+H\left(2+\frac{\rho_{\mathrm{m}}}{\rho_{\mathrm{m}}+\frac{2}{3}q^{2}}\right)\dot{\mathcal{R}}_{\delta\rho_{\mathrm{m}}}-\frac{1}{3}q^{2}\left(\frac{\rho_{\mathrm{m}}}{\rho_{\mathrm{m}}+\frac{2}{3}q^{2}}\right)\mathcal{R}_{\delta\rho_{\mathrm{m}}}=0$
(12)
$\displaystyle\ddot{\mathcal{R}}_{\delta\rho_{\mathrm{r}}}+H\dot{\mathcal{R}}_{\delta\rho_{\mathrm{r}}}+\frac{1}{3}q^{2}\mathcal{R}_{\delta\rho_{\mathrm{r}}}=F$
(13)

with

$F=-2\left(\frac{\rho_{\mathrm{m}}}{\rho_{\mathrm{m}}+\frac{2}{3}q^{2}}\right)\left(H\dot{\mathcal{R}}_{\delta\rho_{\mathrm{m}}}-\frac{1}{3}q^{2}\mathcal{R}_{\delta\rho_{\mathrm{m}}}\right).$
(14)

With the adiabatic condition

$\displaystyle\mathcal{R}_{\delta\rho_{\mathrm{m}}}$
$\displaystyle=\mathcal{R}_{\delta\rho_{\mathrm{r}}},$ (15)
$\displaystyle{\dot{\mathcal{R}}}_{\delta\rho_{\rm m}}$
$\displaystyle={\dot{\mathcal{R}}}_{\delta\rho_{\rm r}}=0$ (16)

we solve Eqs. (12) and (13) analytically.

The phase transition of thermal inflation is controlled by the temperature of
the radiation of the universe. Hence, the curvature perturbation at the end of
thermal inflation is equal to $\mathcal{R}_{\delta\rho_{\mathrm{r}}}$ and
gives the power spectrum by

$\mathinner{P\mathopen{\left(k\right)}}=\mathinner{P_{\rm
pri}\mathopen{\left(k\right)}}\mathinner{\mathcal{T}^{2}\mathopen{\left(\frac{k}{k_{\rm
b}}\right)}}$ (17)

where $\mathinner{P_{\rm pri}\mathopen{\left(k\right)}}$ is the power spectrum
of the primordial inflation and

$\mathinner{\mathcal{T}\mathopen{\left(\frac{k}{k_{\rm
b}}\right)}}=\cos\left[\left(\frac{k}{k_{\rm
b}}\right)\int_{0}^{\infty}\frac{\mathinner{d\alpha}}{\sqrt{\alpha(2+\alpha^{3})}}\right]\\\
+6\left(\frac{k}{k_{\rm
b}}\right)\int_{0}^{\infty}\frac{\mathinner{d\gamma}}{\gamma^{3}}\int_{0}^{\gamma}\mathinner{d\beta}\left(\frac{\beta}{2+\beta^{3}}\right)^{3/2}\sin\left[\left(\frac{k}{k_{\rm
b}}\right)\int_{\gamma}^{\infty}\frac{\mathinner{d\alpha}}{\sqrt{\alpha(2+\alpha^{3})}}\right]~{}.$
(18)

is the transfer function summarizing the effects of thermal inflation eras on
the evolution of perturbations.

In this paper, we extend the analysis of [23] to $t\lesssim t_{\rm a}$ and
$k\gtrsim k_{\rm a}$ (in contrast to Eqs. (2) and (3)) by treating the moduli
as simple matter, $\rho_{m}\propto a^{-3}$. We study Eqs. (II) and (II) in the
limit $k\gg k_{\rm a}$ in section III and IV. We find numerical solutions for
Eqs. (II) and (II) in section V (in comparison with the analytic solution of
Eq. (18) for Eqs. (12) and (13) in [23]).

##  III Perturbation for thermal inflation plus radiation

To model thermal inflation, we first consider the minimal system of radiation
and vacuum

$\rho=\rho_{\rm r}+V_{0}\,.$ (19)

The curvature perturbation on uniform radiation density hypersurfaces
satisfies

$\displaystyle{\ddot{\mathcal{R}}}_{\delta{\rho}_{\mathrm{r}}}+H\left(1+\frac{4\rho_{\rm
r}}{2\rho_{\rm
r}+q^{2}}\right){\dot{\mathcal{R}}}_{\delta{\rho}_{\mathrm{r}}}+\frac{1}{3}q^{2}\left(1-\frac{4\rho_{\rm
r}}{2\rho_{\rm r}+q^{2}}\right)\mathcal{R}_{\delta\rho_{\mathrm{r}}}\simeq 0$
(20)

from Eq. (II) and we take a growing mode initial condition

$\mathcal{R}_{\delta\rho_{\mathrm{r}}}=\mathcal{R}_{0}{\quad,\quad}{\dot{\mathcal{R}}}_{\delta\rho_{\rm
r}}=0$ (21)

well outside the horizon during radiation domination.

We define the moment when thermal inflation begins, $t_{\rm b}$, by

$\mathinner{\ddot{a}\mathopen{\left(t_{\rm b}\right)}}\equiv 0~{},$ (22)

or,

$\mathinner{\rho_{\rm r}\mathopen{\left(t_{\rm b}\right)}}=V_{0}~{}.$ (23)

From the characteristic scale at $t=t_{\rm b}$

$k_{\rm b}=a_{\rm b}H_{\rm b}$ (24)

we introduce new parameters

$\beta\equiv\frac{a}{a_{\rm b}}$ (25) $\lambda\equiv\frac{k}{k_{\rm b}}~{},$
(26)

the energy density is parametrized by

$\mathinner{\rho\mathopen{\left(\beta\right)}}=V_{0}\left(1+\beta^{-4}\right)~{},$
(27)

and Eq. (20) becomes

$\beta^{2}\frac{d^{2}\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\beta^{2}}+2\left(1-\frac{1}{1+\beta^{4}}+\frac{1}{1+\frac{\lambda^{2}}{3}\beta^{2}}\right)\beta\frac{d\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\beta}+\frac{2\lambda^{2}\beta^{2}}{3\left(1+\beta^{4}\right)}\left(1-\frac{2}{1+\frac{\lambda^{2}}{3}\beta^{2}}\right)\mathcal{R}_{\delta\rho_{\mathrm{r}}}=0~{}.$
(28)

Now we study Eq. (28) on large scale $\lambda\ll 1$ ($k\ll k_{\rm b}$) and on
small scales $\lambda\gg 1$ ($k\gg k_{\rm b}$) as follows.

###  III.1 Large scales $\lambda\ll 1$

For modes that remain outside the horizon, Eq. (28) reduces to

$\beta^{2}\frac{d^{2}\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\beta^{2}}+2\left(2-\frac{1}{1+\beta^{4}}\right)\beta\frac{d\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\beta}\simeq\frac{2\lambda^{2}\beta^{2}}{3}\left(\beta\frac{d\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\beta}+\frac{1}{1+\beta^{4}}\mathcal{R}_{\delta\rho_{\mathrm{r}}}\right)~{},$
(29)

which has a solution

$\mathcal{R}_{\delta\rho_{\mathrm{r}}}=\mathcal{R}_{0}\left[1+\frac{2}{3}\lambda^{2}\int_{0}^{\beta}\frac{\xi^{2}d\xi}{\sqrt{1+\xi^{4}}}\int_{\xi}^{\beta}\frac{d\eta}{\eta^{2}\sqrt{1+\eta^{4}}}+\mathinner{\mathcal{O}\mathopen{\left(\lambda^{4}\right)}}\right]~{}.$
(30)

###  III.2 Small scales $\lambda\gg 1$

For modes that enter the horizon during radiation domination, we solve Eq.
(28) in two overlapping regimes:

Radiation domination

For $\beta\ll 1$, the modes start well outside the horizon during radiation
domination and end well inside the horizon during radiation domination.Eq.
(28) reduces to

$\beta^{2}\frac{d^{2}\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\beta^{2}}+\left(\frac{2}{1+\frac{1}{3}\lambda^{2}\beta^{2}}\right)\beta\frac{d\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\beta}+\frac{2}{3}\lambda^{2}\beta^{2}\left(1-\frac{2}{1+\frac{1}{3}\lambda^{2}\beta^{2}}\right)\mathcal{R}_{\delta\rho_{\mathrm{r}}}=0~{},$
(31)

which has a solution

$\mathcal{R}_{\delta\rho_{\mathrm{r}}}=A_{1}\left[\frac{\sqrt{6}}{\lambda\beta}\sin\left(\sqrt{\frac{2}{3}}\lambda\beta\right)-\cos\left(\sqrt{\frac{2}{3}}\lambda\beta\right)\right]\\\
+B_{1}\left[\frac{\sqrt{6}}{\lambda\beta}\cos\left(\sqrt{\frac{2}{3}}\lambda\beta\right)+\sin\left(\sqrt{\frac{2}{3}}\lambda\beta\right)\right]~{}.$
(32)

Matching to the initial condition of Eq. (21) at $\beta\ll\lambda^{-1}$ gives

$A_{1}=\mathcal{R}_{0}{\quad,\quad}B_{1}=0~{}.$ (33)

Well inside the horizon during radiation domination to vacuum domination

For $\beta\gg\lambda^{-1}$, the modes start well inside the horizon during
radiation domination and end well outside the horizon during vacuum
domination. Eq. (28) reduces to

$\left(1+\beta^{4}\right)\frac{d^{2}\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\beta^{2}}+2\beta^{3}\frac{d\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\beta}+\frac{2\lambda^{2}}{3}\mathcal{R}_{\delta\rho_{\mathrm{r}}}=0~{},$
(34)

which has a solution

$\displaystyle\mathcal{R}_{\delta\rho_{\mathrm{r}}}=A_{2}\cos\left[\lambda\sqrt{\frac{2}{3}}\int_{0}^{\beta}\frac{1}{\sqrt{1+\xi^{4}}}d\xi\right]+B_{2}\sin\left[\lambda\sqrt{\frac{2}{3}}\int_{0}^{\beta}\frac{1}{\sqrt{1+\xi^{4}}}d\xi\right]~{}.$
(35)

Matching to Eq. (32) at $\lambda^{-1}\ll\beta\ll 1$ gives

$A_{2}=-A_{1}=-\mathcal{R}_{0}{\quad,\quad}B_{2}=B_{1}=0~{}.$ (36)

##  IV Perturbations for thermal inflation plus radiation and matter

To model the primordial thermal bath, the moduli generated thereby, and
thermal inflation, we consider a system of radiation, matter and vacuum

$\rho=\rho_{\rm r}+\rho_{\rm m}+V_{0}=\rho_{0}\left(\frac{a_{\rm
a}}{a}\right)^{4}+\rho_{0}\left(\frac{a_{\rm a}}{a}\right)^{3}+V_{0}~{},$ (37)

where subscript “${\rm a}$” indicates the moment when the matter energy
density surpasses the radiation energy density. We describe the changes of
three phases in the thermal inflation scenario by Eq. (37). The primordial
inflation is followed by radiation-dominated phase, and then moduli dominates
over the universe leaving a matter phase. Thermal inflation occurs to dilute
the moduli matter and drives a vacuum domination. However, Eq. (37) could lead
to two phases — radiation and vacuum phases for $V_{0}/\rho_{o}\gg 1$ as
follows.

Figure 2: The ratios $k_{\rm a}/k_{\rm b}$ and $a_{\rm a}/a_{\rm b}$ as a
function of $\upsilon$.

We introduce a useful parameter

$\upsilon\equiv\frac{V_{0}}{\rho_{0}}~{}$ (38)

which is the ratio of vacuum energy to radiation at matter-radiation equality.
The moment when thermal inflation begins, $t_{\rm b}$, is found by

$\mathinner{\ddot{a}\mathopen{\left(t_{\rm b}\right)}}\equiv 0~{},$ (39)

giving

$2\left(\frac{a_{\rm a}}{a_{\rm b}}\right)^{4}+\left(\frac{a_{\rm a}}{a_{\rm
b}}\right)^{3}-2\upsilon=0$ (40)

and

$\left(2+\upsilon\right)\left(\frac{k_{\rm b}}{k_{\rm
a}}\right)^{2}=\left(\frac{a_{\rm b}}{a_{\rm
a}}\right)^{2}\left[\left(\frac{a_{\rm a}}{a_{\rm
b}}\right)^{4}+\left(\frac{a_{\rm a}}{a_{\rm
b}}\right)^{3}+\upsilon\right]~{},$ (41)

where

$k_{\rm a}=a_{\rm a}H_{\rm a}$ (42)

and $k_{\rm b}=a_{\rm b}H_{\rm b}$ is already given in Eq. (24). From Eqs.
(40) and (41), the relations between $\upsilon$, $a_{\rm a}/a_{\rm b}$ and
$k_{\rm a}/k_{\rm b}$ are fixed by the value of $\upsilon$ in FIG. 2.

Figure 3:  A diagram showing the evolution of horizon scales for $\upsilon\ll
1$. The characteristic scales $k_{\rm a}$ and $k_{\rm b}$ are shown with the
phase evolution.

For $\upsilon\ll 1$, Eqs. (40) and (41) produces asymptotic formulae

$\frac{a_{\rm a}}{a_{\rm
b}}\simeq\left(2\upsilon\right)^{\frac{1}{3}}{\quad,\quad}\frac{k_{\rm
a}}{k_{\rm
b}}\simeq\left(\frac{2^{5}}{3^{3}\upsilon}\right)^{\frac{1}{6}}~{},$ (43)

implying $a_{\rm a}\ll a_{\rm b}$ and $k_{\rm b}\ll k_{\rm a}$. In this case,
we have a radiation phase initially at $a\ll a_{\rm a}$, a matter phase at
$a_{\rm a}<a<a_{\rm b}$, and then a vacuum phase at $a>a_{\rm b}$ in FIG. 3.
These three phases mimic the eras after primordial inflation in thermal
inflation scenario. In this limit, Eq. (37) gives $\mathinner{\rho_{\rm
r}\mathopen{\left(a_{\rm b}\right)}}\ll\mathinner{\rho_{\rm
m}\mathopen{\left(a_{\rm b}\right)}}$, and it reduces to the case of [23].

Figure 4: A diagram showing the evolution of horizon scales for $\upsilon\gg
1$. The characteristic scales $k_{\rm a}$ and $k_{\rm b}$ are shown with the
phase evolution. Note that matter-domination phase from Eq. (37) does not
exist.

For $\upsilon\gg 1$, Eqs. (40) and (41) produces asymptotic relations

$\frac{a_{\rm a}}{a_{\rm
b}}\simeq\upsilon^{\frac{1}{4}}{\quad,\quad}\frac{k_{\rm a}}{k_{\rm
b}}\simeq\left(\frac{\upsilon}{4}\right)^{\frac{1}{4}}~{},$ (44)

implying $a_{\rm b}\ll a_{\rm a}$ and $k_{\rm b}\ll k_{\rm a}$. In this case,
we have a radiation phase initially at $a\ll a_{\rm b}$ which turns into a
vacuum phase at $a=a_{\rm b}$. At $a=a_{\rm a}$, however, Eq. (37) gives
$\rho_{\rm m}+\rho_{\rm r}\ll V_{0}$ leaving a vacuum phase. Hence, Eq. (37)
describes only two phases of radiation and vacuum in FIG. 4. In this limit,
Eq. (37) gives $\mathinner{\rho_{\rm r}\mathopen{\left(a_{\rm
b}\right)}}\gg\mathinner{\rho_{\rm m}\mathopen{\left(a_{\rm b}\right)}}$, and
it clearly becomes the case of thermal inflation plus radiation discussed in
the previous section.

We take the adiabatic growing mode as the initial condition

$\displaystyle\mathcal{R}_{\delta\rho_{\mathrm{m}}}$
$\displaystyle=\mathcal{R}_{\delta\rho_{\mathrm{r}}}=\mathcal{R}_{0},$ (45)
$\displaystyle{\dot{\mathcal{R}}}_{\delta\rho_{\rm m}}$
$\displaystyle={\dot{\mathcal{R}}}_{\delta\rho_{\rm r}}=0$ (46)

well outside the horizon during radiation domination.

We rewrite the equations by introducing new parameters

$\displaystyle\alpha$ $\displaystyle\equiv\frac{a}{a_{\rm a}}~{},$ (47)
$\displaystyle\kappa$ $\displaystyle\equiv\frac{k}{k_{\rm a}}~{}.$ (48)

Then Eq. (II) becomes

$\alpha^{2}\frac{d^{2}\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\alpha^{2}}+2\left(1+\frac{1}{1+\frac{3}{4}\alpha+\frac{2+\upsilon}{6}\kappa^{2}\alpha^{2}}-\frac{1+\frac{3}{4}\alpha}{1+\alpha+\upsilon\alpha^{4}}\right)\alpha\frac{d\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\alpha}\\\
+\left(\frac{2+\upsilon}{3}\kappa^{2}\alpha^{2}\right)\left(\frac{1}{1+\alpha+\upsilon\alpha^{4}}\right)\left(1-\frac{2}{1+\frac{3}{4}\alpha+\frac{2+\upsilon}{6}\kappa^{2}\alpha^{2}}\right)\mathcal{R}_{\delta\rho_{\mathrm{r}}}\\\
=-\frac{3}{2}\left(\frac{\alpha}{1+\frac{3}{4}\alpha+\frac{2+\upsilon}{6}\kappa^{2}\alpha^{2}}\right)\left[\alpha\frac{d\mathcal{R}_{\delta\rho_{\mathrm{m}}}}{d\alpha}-\left(\frac{2+\upsilon}{3}\kappa^{2}\alpha^{2}\right)\left(\frac{1}{1+\alpha+\upsilon\alpha^{4}}\right)\mathcal{R}_{\delta\rho_{\mathrm{m}}}\right]~{},$
(49)

and Eq. (II) becomes

$\alpha^{2}\frac{d^{2}\mathcal{R}_{\delta\rho_{\mathrm{m}}}}{d\alpha^{2}}+\left(3+\frac{\frac{3}{4}\alpha}{1+\frac{3}{4}\alpha+\frac{2+\upsilon}{6}\kappa^{2}\alpha^{2}}-\frac{2+\frac{3}{2}\alpha}{1+\alpha+\upsilon\alpha^{4}}\right)\alpha\frac{d\mathcal{R}_{\delta\rho_{\mathrm{m}}}}{d\alpha}\\\
-\left(\frac{2+\upsilon}{3}\kappa^{2}\alpha^{2}\right)\left(\frac{1}{1+\alpha+\upsilon\alpha^{4}}\right)\left(\frac{\frac{3}{4}\alpha}{1+\frac{3}{4}\alpha+\frac{2+\upsilon}{6}\kappa^{2}\alpha^{2}}\right)\mathcal{R}_{\delta\rho_{\mathrm{m}}}\\\
=-\left(\frac{1}{1+\frac{3}{4}\alpha+\frac{2+\upsilon}{6}\kappa^{2}\alpha^{2}}\right)\left[\alpha\frac{d\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\alpha}-\left(\frac{2+\upsilon}{3}\kappa^{2}\alpha^{2}\right)\left(\frac{1}{1+\alpha+\upsilon\alpha^{4}}\right)\mathcal{R}_{\delta\rho_{\mathrm{r}}}\right]~{}.$
(50)

Now we study the asymptotic behaviors of Eqs. (49) and (50) for both
$\upsilon\ll 1$ and $\upsilon\gg 1$ on large scale $\kappa\ll 1$ ($k\ll k_{\rm
a}$) and on small scales $\kappa\gg 1$ ($k\gg k_{\rm a}$) as follows.

###  IV.1 Large scales $\kappa\ll 1$

We consider the cases of $\upsilon\ll 1$ and $\upsilon\gg 1$ separately in the
large scale limit. For

$\upsilon\ll 1~{},$ (51)

and for the late time

$\alpha\gg 1~{},$ (52)

Eqs. (II) and (II) become

$\left(\alpha+\upsilon\alpha^{4}\right)\frac{d^{2}\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\alpha^{2}}+\frac{1}{2}\left(1+4\upsilon\alpha^{3}\right)\frac{d\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\alpha}+\frac{2}{3}\kappa^{2}\mathcal{R}_{\delta\rho_{\mathrm{r}}}\\\
=-\frac{3}{2}\left(\frac{1}{\frac{3}{4}\alpha+\frac{1}{3}\kappa^{2}\alpha^{2}}\right)\left[\left(\alpha+\upsilon\alpha^{4}\right)\frac{d\mathcal{R}_{\delta\rho_{\mathrm{m}}}}{d\alpha}-\frac{2}{3}\kappa^{2}\alpha\mathcal{R}_{\delta\rho_{\mathrm{m}}}\right]$
(53)
$\alpha^{2}\frac{d^{2}\mathcal{R}_{\delta\rho_{\mathrm{m}}}}{d\alpha^{2}}+\left(3+\frac{1}{1+\frac{4}{9}\kappa^{2}\alpha}-\frac{\frac{3}{2}}{1+\upsilon\alpha^{3}}\right)\alpha\frac{d\mathcal{R}_{\delta\rho_{\mathrm{m}}}}{d\alpha}\\\
-\left(\frac{\frac{2}{3}\kappa^{2}\alpha}{1+\upsilon\alpha^{3}}\right)\left(\frac{1}{1+\frac{4}{9}\kappa^{2}\alpha}\right)\mathcal{R}_{\delta\rho_{\mathrm{m}}}=0~{},$
(54)

which have solutions

$\mathcal{R}_{\delta\rho_{\mathrm{r}}}=A\cos\left[\sqrt{\frac{2}{3}}\kappa\int_{\alpha_{i}}^{\alpha}\frac{d\xi}{\sqrt{\xi+\upsilon\xi^{4}}}\right]+B\sin\left[\sqrt{\frac{2}{3}}\kappa\int_{\alpha_{i}}^{\alpha}\frac{d\xi}{\sqrt{\xi+\upsilon\xi^{4}}}\right]\\\
+\frac{9}{2}\int_{\alpha_{i}}^{\alpha}d\xi\frac{\sqrt{\frac{2}{3}}\kappa}{\sqrt{\xi+\upsilon\xi^{4}}}\sin\left[\sqrt{\frac{2}{3}}\kappa\int_{\xi}^{\alpha}\frac{d\sigma}{\sqrt{\sigma+\upsilon\sigma^{4}}}\right]\left[\frac{\mathinner{\mathcal{R}_{\delta\rho_{\mathrm{m}}}\mathopen{\left(\xi\right)}}-C}{\kappa^{2}\xi}\right]~{},$
(55)

where $\alpha_{i}\ll 1$ is the initial value, and

$\displaystyle\mathcal{R}_{\delta\rho_{\mathrm{m}}}=C\left[1+\frac{2}{3}\kappa^{2}\sqrt{\frac{1+\upsilon\alpha^{3}}{\alpha^{3}}}\int_{0}^{\alpha}\left(\frac{\xi}{1+\upsilon\xi^{3}}\right)^{\frac{3}{2}}d\xi\right]+D\sqrt{\frac{1+\upsilon\alpha^{3}}{\alpha^{3}}}~{}.$
(56)

Matching to the initial conditions of Eqs. (45) and (46) at
$\alpha\ll\kappa^{-2}$ gives

$A=\mathcal{R}_{0}{\quad,\quad}B=0{\quad,\quad}C=\mathcal{R}_{0}{\quad,\quad}D=0$
(57)

reproducing the result of [23].

For $\upsilon\gg 1$, we return to the case of thermal inflation plus
radiation, and Eqs. (II) and (II) reduce to Eq. (20) giving Eq. (30) as the
proper solution.

###  IV.2 Small scales $\kappa\gg 1$

We study the cases of $\upsilon\ll 1$ and $\upsilon\gg 1$ together in the
small scale limit. For modes that enter the horizon during radiation
domination, we solve Eqs. (49) and (50) in two overlapping regimes:

Radiation domination

For $\alpha\ll\mathinner{{\rm
min}\mathopen{\left(1,\upsilon^{-\frac{1}{4}}\right)}}$, Eq. (49) becomes

$\alpha^{2}\frac{d^{2}\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\alpha^{2}}+\left(\frac{2}{1+\frac{2+\upsilon}{6}\kappa^{2}\alpha^{2}}\right)\alpha\frac{d\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\alpha}+\frac{2+\upsilon}{3}\kappa^{2}\alpha^{2}\left(1-\frac{2}{1+\frac{2+\upsilon}{6}\kappa^{2}\alpha^{2}}\right)\mathcal{R}_{\delta\rho_{\mathrm{r}}}=0~{},$
(58)

whose solution is given by

$\mathcal{R}_{\delta\rho_{\mathrm{r}}}=A_{1}\left[\sqrt{\frac{3}{2+\upsilon}}\left(\frac{2}{\kappa\alpha}\right)\sin\left(\sqrt{\frac{2+\upsilon}{3}}\kappa\alpha\right)-\cos\left(\sqrt{\frac{2+\upsilon}{3}}\kappa\alpha\right)\right]\\\
+B_{1}\left[\sqrt{\frac{3}{2+\upsilon}}\left(\frac{2}{\kappa\alpha}\right)\cos\left(\sqrt{\frac{2+\upsilon}{3}}\kappa\alpha\right)+\sin\left(\sqrt{\frac{2+\upsilon}{3}}\kappa\alpha\right)\right]~{}.$
(59)

Matching to the initial conditions of Eqs. (45) and (46) at
$\alpha\ll\mathinner{{\rm
min}\mathopen{\left(\kappa^{-1},\kappa^{-1}\upsilon^{-\frac{1}{2}}\right)}}$
gives

$A_{1}=\mathcal{R}_{0}{\quad,\quad}B_{1}=0~{}.$ (60)

Well inside the horizon during radiation domination to vacuum domination

For $\alpha\gg\mathinner{{\rm
min}\mathopen{\left(\kappa^{-1},\kappa^{-1}\upsilon^{-\frac{1}{2}}\right)}}$,
Eq. (49) becomes

$\left(1+\alpha+\upsilon\alpha^{4}\right)\frac{d^{2}\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\alpha^{2}}+\frac{1}{2}\left(1+4\upsilon\alpha^{3}\right)\frac{d\mathcal{R}_{\delta\rho_{\mathrm{r}}}}{d\alpha}+\left(\frac{2+\upsilon}{3}\right)\kappa^{2}\mathcal{R}_{\delta\rho_{\mathrm{r}}}=0~{},$
(61)

whose solution is given by

$\mathcal{R}_{\delta\rho_{\mathrm{r}}}=A_{2}\cos\left(\sqrt{\frac{2+\upsilon}{3}}\kappa\int_{0}^{\alpha}\frac{d\xi}{\sqrt{1+\xi+\upsilon\xi^{4}}}\right)\\\
+B_{2}\sin\left(\sqrt{\frac{2+\upsilon}{3}}\kappa\int_{0}^{\alpha}\frac{d\xi}{\sqrt{1+\xi+\upsilon\xi^{4}}}\right)~{}.$
(62)

Matching to Eq. (59) at $\mathinner{{\rm
min}\mathopen{\left(\kappa^{-1},\kappa^{-1}\upsilon^{-\frac{1}{2}}\right)}}\ll\alpha\ll\mathinner{{\rm
min}\mathopen{\left(1,\upsilon^{-\frac{1}{4}}\right)}}$ gives

$A_{2}=-A_{1}=-\mathcal{R}_{0}{\quad,\quad}B_{2}=B_{1}=0~{}.$ (63)

For $\upsilon\gg 1$, we can show that Eq. (62) becomes Eq. (35) by using the
relation of Eq. (43).

In the next section, we construct the asymptotic forms of transfer function
from the primordial inflationary to thermal inflation power spectrum by using
the results of Eqs. (30) and (62) for $\upsilon\gg 1$ and Eqs. (55) and (62)
for $\upsilon\ll 1$ to study the numerical solution of Eqs. (49) and (50) for
various values of $\upsilon$.

##  V Transfer functions

Figure 5: Transfer functions of the thermal inflation scenarios with various
values of $\upsilon$ as a function of $k/k_{\rm b}$.

We take the adiabatic condition to calculate the density perturbation in Eqs.
(45) and (46). The effects of thermal inflation is summarized by the curvature
perturbation on the radiation density hypersurfaces
$\mathcal{R}_{\delta\rho_{\mathrm{r}}}$ as

$\mathinner{\mathcal{T}\mathopen{\left(k\right)}}\equiv\frac{\mathinner{\mathcal{R}_{\delta\rho_{\mathrm{r}}}\mathopen{\left(k,\infty\right)}}}{\mathinner{\mathcal{R}_{\delta\rho_{\mathrm{r}}}\mathopen{\left(k,0\right)}}}~{}.$
(64)

For thermal inflation and radiation, we find the asymptotic form of the
transfer function from Eqs. (30) and (35). On large scale, it is almost scale-
invariant and has slight enhancement at $k\simeq k_{\rm b}$,

$\mathinner{\mathcal{T}\mathopen{\left(k\right)}}\xrightarrow{k\ll k_{\rm
b}}1+\mu_{0}\left(\frac{k}{k_{\rm b}}\right)^{2},$ (65)

where

$\mu_{0}=\int_{0}^{\infty}\frac{\xi^{2}d\xi}{\sqrt{1+\xi^{4}}}\int_{\xi}^{\infty}\frac{d\eta}{\eta^{2}\sqrt{1+\eta^{4}}}\simeq
0.2393~{}.$ (66)

On small scales, it is sinusoidal with an amplitude of unity,

$\mathinner{\mathcal{T}\mathopen{\left(k\right)}}\xrightarrow{k\gg k_{\rm
b}}-\cos\left[\mu_{1}\left(\frac{k}{k_{\rm b}}\right)\right]~{},$ (67)

where

$\mu_{1}=\sqrt{\frac{2}{3}}\int_{0}^{\infty}\frac{d\beta}{\sqrt{1+\beta^{4}}}=\frac{\mathinner{\Gamma\mathopen{\left(\frac{1}{4}\right)}}^{2}}{2\sqrt{6\pi}}\simeq
1.5139~{}.$ (68)

Figure 6: Transfer functions of the thermal inflation scenarios with various
values of $\upsilon\lesssim 1$ as a function of $k/k_{\rm a}$.

For thermal inflation, radiation and matter, the transfer function for
$\upsilon\ll 1$ becomes

$\mathinner{\mathcal{T}\mathopen{\left(k\right)}}=\cos\left[\left(\frac{k}{k_{{\rm
b}}}\right)\int_{0}^{\infty}\frac{d\xi}{\sqrt{\xi(2+\xi^{3})}}\right]\\\
+6\left(\frac{k}{k_{{\rm
b}}}\right)\int_{0}^{\infty}\frac{d\eta}{\eta^{3}}\int_{0}^{\eta}d\zeta\left(\frac{\zeta}{2+\zeta^{3}}\right)^{\frac{3}{2}}\sin\left[\left(\frac{k}{k_{{\rm
b}}}\right)\int_{\eta}^{\infty}\frac{d\xi}{\sqrt{\xi(2+\xi^{3})}}\right]$ (69)

from Eq. (55), which is exactly the same form of [23].

On large scales that remain outside the horizon, the transfer function goes to

$\mathinner{\mathcal{T}\mathopen{\left(k\right)}}\longrightarrow
1+\nu_{0}\left(\frac{k}{k_{\rm b}}\right)^{2}~{},$ (70)

where

$\nu_{0}=\int_{0}^{\infty}\mathinner{d\alpha}\left(\frac{\alpha}{2+\alpha^{3}}\right)^{\frac{3}{2}}=\frac{2^{7/3}\pi^{3/2}}{3^{3/2}\mathinner{\Gamma\mathopen{\left(\frac{1}{6}\right)}}\mathinner{\Gamma\mathopen{\left(\frac{1}{3}\right)}}}\simeq
0.3622~{}.$ (71)

On smaller scales that enter the horizon during matter domination, $k_{\rm
b}\ll k\ll k_{\rm a}$, the transfer function goes to

$\mathinner{\mathcal{T}\mathopen{\left(k\right)}}\longrightarrow-\frac{1}{5}\cos\left[\nu_{1}\left(\frac{k}{k_{\rm
b}}\right)\right]~{},$ (72)

where

$\nu_{1}=\int_{0}^{\infty}\frac{\mathinner{d\alpha}}{\sqrt{\alpha(2+\alpha^{3})}}=\frac{\mathinner{\Gamma\mathopen{\left(\frac{1}{6}\right)}}\mathinner{\Gamma\mathopen{\left(\frac{1}{3}\right)}}}{2^{1/3}3\sqrt{\pi}}\simeq
2.2258~{}.$ (73)

On much smaller scales that enter the horizon during radiation domination,
$k\gg k_{\rm a}$, Eq. (62) gives the transfer function as

$\mathinner{\mathcal{T}\mathopen{\left(k\right)}}\longrightarrow-\cos\left[\mathinner{\tau_{0}\mathopen{\left(\upsilon\right)}}\sqrt{\frac{2+\upsilon}{3}}\left(\frac{k}{k_{\rm
a}}\right)\right]\simeq\begin{cases}-\cos\left[\nu_{1}\left(\frac{k}{k_{\rm
b}}\right)\right]&\text{for $\upsilon\ll 1$}\\\
-\cos\left[\mu_{1}\left(\frac{k}{k_{\rm b}}\right)\right]&\text{for
$\upsilon\gg 1$}\end{cases}~{},$ (74)

where

$\mathinner{\tau_{0}\mathopen{\left(\upsilon\right)}}=\int_{0}^{\infty}\frac{d\alpha}{\sqrt{1+\alpha+\upsilon\alpha^{4}}}\simeq\begin{cases}2^{\frac{1}{3}}\nu_{1}\upsilon^{-\frac{1}{6}}-2+\mathinner{\mathcal{O}\mathopen{\left(\upsilon^{\frac{1}{6}}\right)}}&\text{for
$\upsilon\ll 1$}\\\
\sqrt{\frac{3}{2}}\mu_{1}\upsilon^{-\frac{1}{4}}-\frac{1}{4}\upsilon^{-\frac{1}{2}}+\mathinner{\mathcal{O}\mathopen{\left(-\upsilon^{\frac{3}{4}}\right)}}&\text{for
$\upsilon\gg 1$}\end{cases}~{}.$ (75)

The transfer functions are plotted with $k/k_{\rm b}$ in FIG. 5, and they are
consistent with the asymptotic forms of Eqs. (65), (70), (72) and (74). In
FIG. 6, the transfer functions for $\upsilon\lesssim 1$ are plotted with
$k/k_{\rm a}$, and we can find the envelop of the amplitude of transfer
functions

$\mathinner{\mathcal{T}_{\rm
env}\mathopen{\left(k\right)}}\simeq\frac{0.2+0.396\left(\frac{k}{k_{\rm
a}}\right)^{1.093}}{1+0.396\left(\frac{k}{k_{\rm
a}}\right)^{1.093}}\longrightarrow\begin{cases}0.2,&k_{\rm b}<k\ll k_{\rm
a}\\\ 0.427,&k\simeq k_{\rm a}\\\ 1,&k_{\rm a}\ll k\end{cases}~{}.$ (76)

Note that the oscillating part with the amplitude of $0.2$ can be found only
for $\upsilon\ll 1$. At $k\simeq k_{\rm a}$, the envelop is more precisely
approximated by

$\mathinner{\mathcal{T}_{\rm env}\mathopen{\left(k\right)}}\longrightarrow
0.427+0.422\log_{10}\left(\frac{k}{k_{\rm a}}\right)~{}.$ (77)

##  VI Conclusion

In this paper, we model density perturbations for thermal inflation by
considering a multi-component system of radiation, matter and vacuum energy
$\rho=\rho_{\rm r}+\rho_{\rm m}+V_{0}$. The transfer function converting the
primordial power spectrum to the model prediction is calculated by tracing the
linear evolution of their curvature perturbations in FIG. 5 and FIG. 6. The
value of $\upsilon\equiv V_{0}/\rho_{0}$ is a key parameter governing the
shape of transfer functions. $\upsilon$ adjusts the relative ratio between
matter and radiation energy densities at the beginning of the second (thermal)
inflation. The system equivalently reduces to the case of [23] for
$\upsilon\rightarrow 0$ and to the case of vacuum and radiation for
$\upsilon\rightarrow\infty$.

For $\upsilon<1.5$ (see FIG. 2), we considered the following three kinds of
modes: first, the largest modes with $k<k_{\rm b}$ are always beyond the
horizon and their curvature perturbation also remains constant leaving
$\mathcal{T}\simeq 1$. Second, the intermediate modes with $k_{\rm b}<k<k_{\rm
a}$ come into the horizon during moduli matter domination and exit the horizon
during thermal inflation. For these modes, perturbations are enhanced by
$\mathcal{T}\simeq 1.3622$ at $k\simeq k_{\rm b}$ and suppressed by
$\left|\mathcal{T}\right|\rightarrow\frac{1}{5}$ on smaller scales $k\gg
k_{\rm b}$. Third, the smallest modes with $k>k_{\rm a}$ enter the horizon
during radiation domination and exit during thermal inflation. These modes are
inside the horizon at radiation domination before the moduli domination and so
behave oscillating with the amplitude $\left|\mathcal{T}\right|\rightarrow 1$.

For $\upsilon\gg 1$, the matter component is negligible compared to the
radiation at the beginning of thermal inflation. The transfer function is
constant on large scales $k<k_{\rm b}$, enhanced as $\mathcal{T}\simeq 1.2393$
at $k\simeq k_{\rm b}$, and then have sinusoidal oscillations with an
amplitude of unity at $k\gg k_{\rm b}$.

Note that the characteristic scale $k_{\rm b}$ is estimated by

$k_{\rm b}\simeq 10^{3}\mathinner{{\rm
Mpc}^{-1}}\left(\frac{e^{20}}{e^{N}}\right)\left(\frac{V_{0}^{\frac{1}{4}}}{10^{7}\mathinner{\mathrm{GeV}}}\right)^{\frac{2}{3}}\left(\frac{T_{\rm
d}}{\mathinner{\mathrm{GeV}}}\right)^{\frac{1}{3}}$ (78)

where $N$ is the e-folds during thermal inflation, the vacuum energy for
thermal inflation is $10^{3}\mathinner{\mathrm{GeV}}\ll
V_{0}^{\frac{1}{4}}<10^{11}\mathinner{\mathrm{GeV}}$, and the reheating
temperature for the radiation domination is
$10^{-2}\mathinner{\mathrm{GeV}}\lesssim T_{\rm d}\lesssim
10^{2}\mathinner{\mathrm{GeV}}$ in [23, 25], and $k_{\rm a}/k_{\rm b}\lesssim
10$ if ${10^{-4}}\lesssim\upsilon\lesssim 10^{4}$ in FIG. 2. Hence, the
changes in the transfer functions according to the value of $\upsilon$ could
be explored by small-scale observations including CMB spectral distortions
[24], the substructure of galaxies [25, 27, 28], and the 21-cm hydrogen line
[25, 29]. Our results can be applied to the calculation for the density
perturbations in multiple inflation scenarios [28, 30, 31].

In this work, we assume the moduli is a simple matter approximated by
$\rho_{\rm m}\propto a^{-3}$ as a baseline model of thermal inflation. The
moduli field oscillates around the minimum affecting perturbations during its
generation and domination, but its detailed dynamics could not be captured in
our assumption. We leave full consideration of moduli dynamics and its
observational implications as a future work.

## Acknowledgements

The authors thank Ewan Stewart for his helpful discussion and advice. HZ
thanks Emre Onur Kahya, and the Department of Physics Engineering at Istanbul
Technical University and Korea Astronomy and Space Science Institute for the
hospitality. This work was supported by the DGIST-UGRP grant. SEH was
supported by the project “Understanding Dark Universe Using Large Scale
Structure of the Universe”, funded by the Ministry of Science.

## References

  * Gliner [1966] E. B. Gliner, Sov. Phys. JETP 22, 378 (1966).
  * Gliner [1970] E. B. Gliner, Sov. Phys. Dokl. 15, 559 (1970).
  * Guth [1981] A. H. Guth, Phys. Rev. D 23, 347 (1981).
  * Linde [1982] A. D. Linde, Phys. Lett. B 108, 389 (1982).
  * Albrecht and Steinhardt [1982] A. Albrecht and P. J. Steinhardt, Phys. Rev. Lett. 48, 1220 (1982).
  * Akrami  _et al._ [2020] Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] .
  * Coughlan  _et al._ [1983] G. D. Coughlan, W. Fischler, E. W. Kolb, S. Raby, and G. G. Ross, Phys. Lett. B131, 59 (1983).
  * Banks  _et al._ [1994] T. Banks, D. B. Kaplan, and A. E. Nelson, Phys. Rev. D49, 779 (1994), arXiv:hep-ph/9308292 [hep-ph] .
  * de Carlos  _et al._ [1993] B. de Carlos, J. A. Casas, F. Quevedo, and E. Roulet, Phys. Lett. B318, 447 (1993), arXiv:hep-ph/9308325 [hep-ph] .
  * Randall and Thomas [1995] L. Randall and S. D. Thomas, Nucl. Phys. B449, 229 (1995), arXiv:hep-ph/9407248 [hep-ph] .
  * Lyth and Stewart [1995] D. H. Lyth and E. D. Stewart, Phys. Rev. Lett. 75, 201 (1995), arXiv:hep-ph/9502417 [hep-ph] .
  * Lyth and Stewart [1996] D. H. Lyth and E. D. Stewart, Phys. Rev. D53, 1784 (1996), arXiv:hep-ph/9510204 [hep-ph] .
  * Yamamoto [1985] K. Yamamoto, Phys. Lett. B161, 289 (1985).
  * Yamamoto [1986] K. Yamamoto, Phys. Lett. B168, 341 (1986).
  * Enqvist  _et al._ [1986] K. Enqvist, D. V. Nanopoulos, and M. Quiros, Phys. Lett. 169B, 343 (1986).
  * Bertolami and Ross [1987] O. Bertolami and G. G. Ross, Phys. Lett. B183, 163 (1987).
  * Ellis  _et al._ [1987] J. R. Ellis, K. Enqvist, D. V. Nanopoulos, and K. A. Olive, Phys. Lett. B188, 415 (1987).
  * Ellis  _et al._ [1989] J. R. Ellis, K. Enqvist, D. V. Nanopoulos, and K. A. Olive, Phys. Lett. B225, 313 (1989).
  * Boyle and Steinhardt [2008] L. A. Boyle and P. J. Steinhardt, Phys. Rev. D 77, 063504 (2008), arXiv:astro-ph/0512014 .
  * Giovannini [2009] M. Giovannini, Class. Quant. Grav. 26, 045004 (2009), arXiv:0807.4317 [astro-ph] .
  * Kuroyanagi  _et al._ [2011] S. Kuroyanagi, K. Nakayama, and S. Saito, Phys. Rev. D 84, 123513 (2011), arXiv:1110.4169 [astro-ph.CO] .
  * Chun  _et al._ [2018] E. J. Chun _et al._ , Int. J. Mod. Phys. A 33, 1842005 (2018), arXiv:1711.02865 [hep-ph] .
  * Hong  _et al._ [2015] S. E. Hong, H.-J. Lee, Y. J. Lee, E. D. Stewart, and H. Zoe, JCAP 06, 002, arXiv:1503.08938 [astro-ph.CO] .
  * Cho  _et al._ [2017] K. Cho, S. E. Hong, E. D. Stewart, and H. Zoe, JCAP 08, 002, arXiv:1705.02741 [astro-ph.CO] .
  * Hong  _et al._ [2017] S. E. Hong, H. Zoe, and K. Ahn, Phys. Rev. D 96, 103515 (2017), arXiv:1706.08049 [astro-ph.CO] .
  * Kodama and Sasaki [1984] H. Kodama and M. Sasaki, Prog. Theor. Phys. Suppl. 78, 1 (1984).
  * Leo  _et al._ [2018] M. Leo, C. M. Baugh, B. Li, and S. Pascoli, JCAP 12, 010, arXiv:1807.04980 [astro-ph.CO] .
  * Enqvist  _et al._ [2020] K. Enqvist, T. Sawala, and T. Takahashi, JCAP 10, 053, arXiv:1905.13580 [astro-ph.CO] .
  * Leo  _et al._ [2020] M. Leo, T. Theuns, C. M. Baugh, B. Li, and S. Pascoli, JCAP 04, 004, arXiv:1909.04641 [astro-ph.CO] .
  * Bae  _et al._ [2018] G. Bae, S. Bae, S. Choe, S. H. Lee, J. Lim, and H. Zoe, Phys. Lett. B 782, 117 (2018), arXiv:1712.04583 [astro-ph.CO] .
  * Allahverdi  _et al._ [2020] R. Allahverdi _et al._ 10.21105/astro.2006.16182 (2020), arXiv:2006.16182 [astro-ph.CO] .