DESY 20-135

KANAZAWA-20-06

TUM-HEP-1279-20

Primordial gravitational waves in a minimal
model of particle physics and cosmology

Andreas<EMAIL_ADDRESS>Ken’ichi
<EMAIL_ADDRESS>and Carlos
<EMAIL_ADDRESS>
|  1 Deutsches Elektronen-Synchrotron DESY,
---|---
| Notkestraße 85, D-22607 Hamburg, Germany
|  2 Institute for Theoretical Physics, Kanazawa University,
| Kakuma-machi, Kanazawa, Ishikawa 920-1192, Japan
|  3 Physik-Department T70, Technische Universität München,
| James-Franck Straße 1, D-85748 Garching, Germany

Abstract

> In this paper we analyze the spectrum of the primordial gravitational waves
> (GWs) predicted in the Standard Model*Axion*Seesaw*Higgs portal inflation
> (SMASH) model, which was proposed as a minimal extension of the Standard
> Model that addresses five fundamental problems of particle physics and
> cosmology (inflation, baryon asymmetry, neutrino masses, strong CP problem,
> and dark matter) in one stroke. The SMASH model has a unique prediction for
> the critical temperature of the second order Peccei-Quinn (PQ) phase
> transition $T_{c}\sim 10^{8}\,\mathrm{GeV}$ up to the uncertainty in the
> calculation of the axion dark matter abundance, implying that there is a
> drastic change in the equation of state of the universe at that temperature.
> Such an event is imprinted on the spectrum of GWs originating from the
> primordial tensor fluctuations during inflation and entering the horizon at
> $T\sim T_{c}$, which corresponds to $f\sim 1\,\mathrm{Hz}$, pointing to a
> best frequency range covered by future space-borne GW interferometers. We
> give a precise estimation of the effective relativistic degrees of freedom
> across the PQ phase transition and use it to evaluate the spectrum of GWs
> observed today. It is shown that the future high sensitivity GW experiment –
> ultimate DECIGO – can probe the nontrivial feature resulting from the PQ
> phase transition in this model.

###### Contents

  1. 1 Introduction
  2. 2 The SMASH model
  3. 3 Equation of state for the Peccei-Quinn phase transition
  4. 4 Spectrum of gravitational waves
  5. 5 Discussion and conclusions
  6. A Improved daisy resummation compatible with thermal decoupling

##  1 Introduction

One of the most robust and model-independent predictions of inflationary
cosmology is a stochastic background of primordial gravitational waves (GWs)
[1, 2] which originates from tensor fluctuations. The power spectrum of tensor
fluctuations, when modes with comoving wave number $k$ exit the horizon during
inflation, is proportional to the square of the Hubble expansion rate
$H=\dot{a}/a$ during inflation:

$\Delta_{t}^{2}(k)=\left.\frac{2}{\pi^{2}}\frac{H^{2}}{M_{P}^{2}}\right|_{k=a_{\rm
inf}H_{\rm inf}}\approx 3.4\times 10^{-12}\left(\frac{H_{\rm
inf}}{10^{13}\,\mathrm{GeV}}\right)^{2},$ (1.1)

where $M_{P}\simeq 2.44\times 10^{18}\,\mathrm{GeV}$ is the reduced Planck
mass and $a_{\rm inf}=a(t_{\rm inf})$ is the cosmic scale factor during
inflation. The amplitude of the primordial GWs is quantified in terms of the
fractional contribution per logarithmic frequency interval ($f=k/(2\pi
a_{0})$, with $a_{0}$ the present scale factor), to the energy in the present
universe,

$\Omega_{\rm gw}(f)\equiv\frac{1}{\rho_{\rm crit}}\frac{d\rho_{\rm
gw}(f)}{d\ln f}=\mathcal{T}_{0}(f)\,\Delta^{2}_{t}(f)\,,$ (1.2)

where $\rho_{\rm gw}$ is the energy density of GWs and $\rho_{\rm
crit}=3H_{0}^{2}M_{P}^{2}$ the critical energy density in terms of today’s
(i.e. at time $t=t_{0}$) Hubble expansion rate $H_{0}=H(t_{0})\approx
3.24\times 10^{-18}\,h$ Hz.

Clearly, a measurement of the overall normalization of the primordial GW
spectrum would yield crucial information about the dynamics of inflation, in
particular its energy scale, $E_{\rm inf}=(3H_{\rm inf}^{2}M_{P}^{2})^{1/4}$.
Moreover, a precise measurement of the frequency dependence of $\Omega_{\rm
gw}(f)$ would provide also information about the post-inflationary expansion
history of the universe. The latter is encoded in Eq. (1.2) in the transfer
function $\mathcal{T}_{0}(f)\equiv\mathcal{T}(t_{0},f)$, which accounts for
the evolution of GWs after their modes reenter the horizon after inflation.
Notably, it is sensitive to the equation of state of the post-inflationary
universe at horizon crossing [3, 4, 5]. In fact, for frequencies corresponding
to modes which reenter the horizon during the post-reheating radiation-
dominated epoch, $\frac{a_{\rm inf}H_{\rm inf}}{2\pi a_{0}}>\frac{a_{\rm
rh}H_{\rm rh}}{2\pi a_{0}}\gg f\gg\frac{a_{\rm eq}H_{\rm eq}}{2\pi
a_{0}}\approx 1.6\times 10^{-17}$ Hz, it can be approximated by [5]

$\mathcal{T}_{0}(f)\,h^{2}\approx 3.2\times 10^{-6}\left(\frac{g_{\ast
s,0}}{3.91}\right)^{\frac{4}{3}}\left(\frac{g_{*\rho}(T_{\rm
hc}(f))}{g_{*s}(T_{\rm hc}(f))}\right)^{\frac{4}{3}}\left[g_{*\rho}(T_{\rm
hc}(f))\right]^{-\frac{1}{3}}\,,$ (1.3)

where $g_{*\rho}(T_{\rm hc}(f))$ and $g_{*s}(T_{\rm hc}(f))$ denote the
effective number of degrees of freedom of the energy and the entropy density
in the hot plasma at the temperature $T_{\rm hc}(f)$ of horizon crossing,
respectively, and $g_{\ast s,0}$ is the value of $g_{\ast s}$ after the
neutrino decoupling. The horizon crossing temperature $T_{\rm hc}$ can be
related to a frequency determined by the Hubble rate at this temperature via

$f=\frac{a_{\rm hc}H_{\rm hc}}{2\pi a_{0}}\approx
1.2\,\mathrm{Hz}\left(\frac{g_{\ast
s,0}}{3.91}\right)^{\frac{1}{3}}\left[\frac{g_{*\rho}(T_{\rm
hc})}{g_{*s}(T_{\rm hc})}\right]^{\frac{1}{2}}\left[{g_{*s}(T_{\rm
hc})}\right]^{\frac{1}{6}}\left(\frac{T_{\rm
hc}}{10^{8}\,\mathrm{GeV}}\right)\,.$ (1.4)

According to this estimate, the spectrum of primordial GWs from inflation is
expected to be almost flat for a huge range of frequencies, $10^{-17}\,{\rm
Hz}\ll f\ll 10^{2}\,{\rm Hz}\,(T_{\rm rh}/10^{10}\,{\rm GeV})$, and of order

$\Omega_{\rm gw}\,h^{2}\approx 1.1\times 10^{-17}\left(\frac{g_{\ast
s,0}}{3.91}\right)^{\frac{4}{3}}\left(\frac{g_{*\rho}(T_{\rm
hc}(f))}{g_{*s}(T_{\rm hc}(f))}\right)^{\frac{4}{3}}\left[g_{*\rho}(T_{\rm
hc}(f))\right]^{-\frac{1}{3}}\left(\frac{H_{\rm
inf}}{10^{13}\,\mathrm{GeV}}\right)^{2}.$ (1.5)

It has features, such as dips and steps, at frequencies corresponding to
temperatures at which the equation of state changes considerably. Notably,
$e^{+}e^{-}$ annihilation, the QCD and the electroweak crossover, occurring at
temperatures around an MeV, 100 MeV and 100 GeV, respectively, are predicted
to show their imprints in the spectrum of primordial GWs at frequencies around
$10^{-11}$ Hz, $10^{-9}$ Hz, and $10^{-6}$ Hz, respectively.

Unfortunately, there are no GW detectors foreseen which would be sensitive
enough to detect the primordial GWs from inflation in this frequency range.
Furthermore, the white dwarf (WD) confusion noise is an obstacle for detecting
primordial GWs at frequencies $f\lesssim 0.1$ Hz. So far, the only possible
instruments to detect directly the primordial gravitational background from
standard inflationary models appear to be future space-borne GW
interferometers, like the Big Bang Observer (BBO) [6, 7, 8, 9] or the Deci-
hertz Interferometer Gravitational Wave Observatory (DECIGO) [10, 11]. In
fact, it has been argued that the ultimate sensitivity of DECIGO may become
comparable to $\Omega_{\rm gw}h^{2}\sim 10^{-20}$[10] at frequencies around 1
Hz – allowing not only for a large statistics detection of the primordial GWs
from inflation, but also for a detailed investigation of possible features
arising from changes in the equation of state at temperatures of the order of
$10^{8}$ GeV.

A particularly well-motivated extension of the Standard Model (SM) which
predicts a sizable change in the equation of state at such temperatures was
developed in Refs. [12, 13, 14]. This minimal extension of the SM – dubbed
Standard Model*Axion*Seesaw*Higgs portal inflation (SMASH) – addresses five
fundamental problems of particle physics and cosmology (the origin of neutrino
masses, the strong CP problem, the nature of dark matter, the generation of
the matter-antimatter asymmetry in the universe, and the nature of the
inflation) in one stroke. Because of its constrained framework, the model
provides definite predictions for various cosmological observables that can be
probed by upcoming axion dark matter and cosmic microwave background (CMB)
experiments. In particular, it predicts a lower bound on the Hubble expansion
rate during inflation, $H_{\rm inf}\gtrsim 10^{13}$ GeV, and a lower bound on
the ratio of the power in tensor to scalar fluctuations, $r\gtrsim 0.004$,
which can be probed by next generation CMB experiments sensitive to the
predicted primordial $B$-mode polarization patterns, such as CMB-S4 [15],
LiteBIRD [16], and the Simons Observatory [17]. Furthermore, it predicts a
reheating temperature around $10^{10}$ GeV, below which the universe has the
equation of state of a thermal relativistic plasma. At around $T_{c}\sim
10^{8}$ GeV, a global Peccei-Quinn (PQ) symmetry gets broken by a second order
phase transition below which a number of particles get massive, leading to a
sizable change in the effective number of relativistic degrees of freedom,
$g_{*\rho}$ and $g_{*s}$, which will show up in a break in the spectrum of the
primordial GWs at around $1\,$Hz.

It is the purpose of this paper to provide a precise prediction of the
primordial GW spectrum expected in SMASH and to confront it with the projected
sensitivity of ultimate DECIGO. This is done as follows. In Sec. 2, we give a
brief review of SMASH, focusing on its inflationary dynamics, its predictions
for the primordial tensor modes, and the nature and the critical temperature
of the PQ phase transition. Section 3 deals with the precise calculation of
the equation of state in SMASH. In Sec. 4, we use the results from the
previous sections to derive the primordial GW spectrum and compare it with the
projected experimental sensitivities. We discuss our results and conclude in
Sec. 5.

##  2 The SMASH model

In the SMASH model [12, 13, 14], the SM is extended by adding a new complex
scalar field $\sigma$ (the PQ field), three singlet neutrinos $N_{i}$, with
$i=1,2,3$, and a vector-like quark $Q$, all charged under a global $U(1)_{\rm
PQ}$ symmetry. The scalar potential in SMASH has the general form

$V(H,\sigma)=\lambda_{H}\left(H^{\dagger}H-\frac{v^{2}}{2}\right)^{2}+\lambda_{\sigma}\left(|\sigma|^{2}-\frac{v_{\sigma}^{2}}{2}\right)^{2}+2\lambda_{H\sigma}\left(H^{\dagger}H-\frac{v^{2}}{2}\right)\left(|\sigma|^{2}-\frac{v_{\sigma}^{2}}{2}\right)\,,$
(2.1)

with $H$ the Higgs doublet, and $\lambda_{H},\lambda_{\sigma}>0$,
$\lambda_{H\sigma}^{2}<\lambda_{H}\lambda_{\sigma}$, such that the electroweak
and PQ symmetry are broken by the vacuum expectation values (VEVs)

$\langle H^{\dagger}H\rangle=v^{2}/2,\hskip
25.83325pt\langle|\sigma|^{2}\rangle=v_{\sigma}^{2}/2\,,$ (2.2)

where $v_{\sigma}\gg v=246$ GeV. The global PQ charges of the fermions and the
hypercharge of the vector quark $Q$ can be chosen such that the only allowed
interactions of the exotic fermions $N_{i},Q$ are

Here we used a four-component notation, with the $N_{i}$ represented by
Majorana spinors, and $L_{i},D_{i}$ denoting the Dirac spinors associated with
the leptons and down quarks of the $i$th generation. The axion arises as a
Goldstone boson associated with the spontaneous breaking of the PQ symmetry
[18, 19, 20]. It can be the main constituent of dark matter if its decay
constant $f_{a}\sim 10^{11}\,\mathrm{GeV}$ [21, 22, 23], where in SMASH we
have $f_{a}=v_{\sigma}$. This new energy scale also provides large Majorana
masses for heavy neutrinos, which explains the smallness of the masses of the
active neutrinos through the seesaw mechanism [24, 25, 26, 27] and leads to
the generation of baryon asymmetry of the universe via the thermal
leptogenesis scenario [28].

In this framework, inflation can arise from the dynamics of the Higgs and PQ
field in the presence of non-minimal couplings to the Ricci scalar $R$,

$S\supset-\int
d^{4}x\sqrt{-g}\,\left[\frac{M^{2}}{2}+\xi_{H}\,H^{\dagger}H+\xi_{\sigma}\,\sigma^{*}\sigma\right]R\,,$
(2.3)

where the mass scale $M$ is related to the actual reduced Planck mass by

$\displaystyle M^{2}_{P}=M^{2}+\xi_{H}v^{2}+\xi_{\sigma}v^{2}_{\sigma}.$ (2.4)

The non-minimal couplings stretch the scalar potential in the Einstein frame,
making it convex and asymptotically flat at large field values. In order to
avoid problems with perturbative unitarity, it is required that
$1\gtrsim\xi_{\sigma}\gg\xi_{H}\geq 0$. Furthermore, the requirement of a
viable reheating demands that $\lambda_{H\sigma}<0$, in order that slow-roll
inflation happens along the line
$h/\rho=\sqrt{-\lambda_{H\sigma}/\lambda_{H}}$, where
$\rho=\sqrt{2}\,|\sigma|$ is the modulus of the PQ field and $h$ the neutral
component of the Higgs doublet in the unitary gauge. Inflation can then be
described in the Einstein frame by a single canonically normalized field
$\chi$ with potential

$\tilde{V}(\chi)=\frac{1}{4}\tilde{\lambda}_{\sigma}\rho(\chi)^{4}\left(1+\xi_{\sigma}\frac{\rho(\chi)^{2}}{M_{P}^{2}}\right)^{-2}\,,\quad\tilde{\lambda}_{\sigma}\equiv\lambda_{\sigma}\left(1-\frac{\lambda_{H\sigma}^{2}}{\lambda_{\sigma}\lambda_{H}}\right)\,.$
(2.5)

The field $\chi$ is the solution of
$\Omega^{2}\,d\chi/d\rho\simeq(b\,\Omega^{2}+6\,\xi_{\sigma}^{2}\,\rho^{2}/M_{P}^{2})^{1/2}$,
with $\Omega\simeq 1+\xi_{\sigma}\,\rho^{2}/M_{P}^{2}$ being the Weyl
transformation into the Einstein frame and
$b=1+|\lambda_{H\sigma}/\lambda_{H}|$. Vacuum stability requires a small value
of $|\lambda_{H\sigma}|\lesssim 10^{-6}$ and consequently $b\approx 1$.

Figure 1: 95% C.L. contours for the parameters of the non-minimally coupled
potential (2.5) giving inflation as constrained by Planck 2018 data at the
pivot scale $0.002$ Mpc-1 [30, 31]. Shown are: the value of the effective
inflationary coupling (left), the predicted tensor-to-scalar ratio (middle),
the Hubble rate during inflation, $H_{\rm inf}$ (right, lower band), and
energy scale during inflation, $\tilde{V}^{1/4}$ (right, upper band), as a
function of the non-minimal coupling parameter $\xi_{\sigma}$. The thicker red
line corresponds to the predictions when accounting for the fact that
reheating in SMASH leads to radiation domination immediately after inflation.
The shaded regions for $\xi_{\sigma}>1$ indicate, approximately, the region
where the predictivity of inflation is threatened by the breakdown of
perturbative unitarity.

The primordial scalar and tensor power spectra, $\Delta^{2}_{s}(k)$ and
$\Delta^{2}_{t}(k)$, can be computed in the slow-roll approximation from the
(potential) slow-roll parameters,

$\epsilon=M_{P}^{2}(\tilde{V}^{\prime}/\tilde{V})\,,\quad\eta=M_{P}^{2}\tilde{V}^{\prime\prime}/\tilde{V}\,,\quad\zeta=M_{P}^{4}\tilde{V}^{\prime}\tilde{V}^{\prime\prime\prime}/\tilde{V}^{2}\,,$
(2.6)

where the primes denote derivatives with respect to $\chi$. Parametrizing the
spectra as

$\displaystyle\Delta^{2}_{s}(k)=A_{s}\left(k/k_{*}\right)^{n_{s}-1+{1/2\,\alpha}\log(k/k_{*})+\cdots}\,,\quad\Delta^{2}_{t}(k)=A_{t}\left(k/k_{*}\right)^{n_{t}+\cdots}\,,$
(2.7)

where all parameters are evaluated at some fiducial scale $k_{*}$, one obtains
at leading order in the slow-roll expansion for the scalar spectral index
$n_{s}$ and its running $\alpha\equiv dn_{s}/d\ln k$,

$n_{s}\simeq
1-6\epsilon+2\eta\,,\quad\alpha\simeq-2\zeta+16\epsilon\,\eta-24\epsilon^{2}\,,$
(2.8)

and for the amplitude of scalar and tensor perturbations, $A_{s}$ and $A_{t}$,

$A_{s}\simeq\frac{1}{24\pi^{2}\epsilon}\frac{\tilde{V}}{M_{P}^{4}}\,,\quad
A_{t}\simeq\frac{2}{3\pi^{2}}\frac{\tilde{V}}{M_{P}^{4}}\,,$ (2.9)

evaluated at the field value corresponding to the time when the scale $k_{*}$
exits the horizon. Therefore, the tensor-to-scalar ratio is given by

$r\equiv\frac{\Delta^{2}_{t}}{\Delta^{2}_{s}}\simeq\frac{A_{t}}{A_{s}}\simeq
16\epsilon\,,$ (2.10)

while the tensor spectral index $n_{t}$ is entirely determined by $r$ through
the consistency relation

$n_{t}\simeq-\frac{r}{8}\,.$ (2.11)

Moreover, the number of e-folds from some initial field value $\rho_{\rm inf}$
until the end of inflation at $\rho_{\rm end}$ can also be approximately
computed analytically:

$N(\chi_{\rm inf},\chi_{\rm end})\simeq\frac{1}{M_{P}}\int_{\chi_{\rm
end}}^{\chi_{\rm
inf}}\frac{d\chi}{\sqrt{2\epsilon}}=\frac{b+6\xi_{\sigma}}{8{M_{P}^{2}}}\left(\rho_{\rm
inf}^{2}-\rho_{\rm
end}^{2}\right)-\frac{3}{4}\ln\left(\frac{M_{P}^{2}+\xi_{\sigma}\rho_{\rm
inf}^{2}}{M_{P}^{2}+\xi_{\sigma}\rho_{\rm end}^{2}}\right)\,.$ (2.12)

Actually, $N(\chi_{\rm inf},\chi_{\rm end})$ can be obtained exactly by
solving the dynamics of the inflaton as a function of the number of e-folds
itself, which is given by equation [29]:

$\displaystyle\frac{d^{2}\chi}{dN^{2}}+3\,\frac{d\chi}{dN}-\frac{1}{2M_{P}^{2}}\left(\frac{d\chi}{dN}\right)^{3}+\left(3M_{P}-\frac{1}{2M_{P}}\left(\frac{d\chi}{dN}\right)^{2}\right)\sqrt{2\epsilon}=0\,,$
(2.13)

and using the condition $\epsilon_{H}\equiv-\dot{{H}}/{H}^{2}=1$ to determine
the value of $\chi$ at the end of inflation.

Figure 1 shows parameters of the non-minimally coupled potential (2.5) giving
inflation as constrained by the CMB observations at Planck 2018 in combination
with lensing, baryon acoustic oscillation (BAO), and BICEP-Keck Array 2015
(BK15) data [30, 31],

$\displaystyle A_{s}$ $\displaystyle=$ $\displaystyle(2.105\pm 0.030)\times
10^{-9}\quad(68\,\%\,\text{C.L., TT,TE,EE+lowE+lensing+BAO}),$ (2.14)
$\displaystyle n_{s}$ $\displaystyle=$ $\displaystyle 0.9665\pm
0.0038\quad(68\,\%\,\text{C.L., TT,TE,EE+lowE+lensing+BAO}),$ (2.15)
$\displaystyle r_{0.002}$ $\displaystyle<$ $\displaystyle
0.058\quad(95\,\%\,\text{C.L., TT,TE,EE+lowE+lensing+BK15+BAO}),$ (2.16)

where the constraints on $A_{s}$ and $n_{s}$ are obtained at the pivot scale
$k_{*}=0.05$ Mpc-1, and that on $r$ is obtained at $k_{*}=0.002$ Mpc-1
(indicated by the subscript). The corresponding tensor power spectrum
predicted in SMASH is shown in Fig. 2, for four different values of the
effective inflationary coupling $\tilde{\lambda}_{\sigma}$. This will be used
later in the determination of the spectrum of primordial GWs.

Figure 2: Tensor power spectrum predicted in SMASH, for four different values
of the effective inflationary coupling $\tilde{\lambda}_{\sigma}$.

Slow-roll inflation ends for $\rho\sim\mathcal{O}(M_{P})$, when the inflaton
field starts undergoing Hubble-damped oscillations in a quartic potential. The
time-averaged stress-energy tensor for these oscillations has the equation of
state of a radiation fluid. Hence, radiation domination starts right after
inflation, and lasts through the phase of reheating in which the oscillating
fields trigger the production of SM particles and the energy of the inflaton
is transferred into the SM plasma. This post-inflationary history in a
radiation-domination era fixes the relation between the scales of the matter
perturbations we observe in the universe today and the size of the primordial
fluctuations which gave rise to them, when they outgrew the Hubble horizon and
became frozen until their later horizon re-entry. This allows to compute the
number of e-folds between horizon crossing and the end of inflation in a
simple way. For a mode with comoving momentum $k$, we have

$N_{e}(k)\simeq\log\frac{a_{\rm eq}\,{H}_{\rm
eq}}{a_{0}\,{H}_{0}}-\frac{1}{4}\log\frac{3\,{H}_{\rm
eq}^{2}}{M_{P}^{2}}-\log\frac{k}{a_{0}\,{H}_{0}}+\frac{1}{2}\log\frac{V_{k}}{M_{P}^{4}}+\frac{1}{4}\log\frac{M_{P}^{4}}{V_{\rm
end}}\,,$ (2.17)

where $V_{k}$ and $V_{\rm end}$ denote, respectively, the energy density at
the time of the mode’s horizon crossing and at the end of inflation. By
equating this expression with the result of integrating (2.13) (or the
simpler, but less accurate expression (2.12)) and using $A_{s}$ to fit
$\tilde{\lambda}_{\sigma}$ we can compute the value of the inflaton field in
the Jordan frame when a given scale $k$ exited the horizon (for a given
$\xi_{\sigma}$) and give a definite prediction for $n_{s}$ and $r$. Using this
expression one gets the thick red lines of Fig. 1.

The Higgs component of the inflaton in SMASH guarantees efficient reheating by
the production of SM gauge bosons. The reheating temperature is predicted to
be around $T_{\rm rh}\sim 10^{10}$ GeV. Such a temperature ensures a thermal
restoration of the PQ symmetry and its spontaneous breakdown later in the
expansion history of the universe. The PQ phase transition in SMASH turns out
to be second order. Its critical temperature can be estimated as [13]

$T_{c}\simeq\frac{2\sqrt{6\lambda_{\sigma}}v_{\sigma}}{\sqrt{8(\lambda_{\sigma}+\lambda_{H\sigma})+\sum_{i}Y_{ii}^{2}+6y^{2}}},$
(2.18)

where $Y_{ii}$ and $y$ are the Yukawa couplings of $N_{i}$ and $Q$,
respectively (see Eq. (2)). Taking into account that the parameters
determining the critical temperature are constrained by the requirements of

  * •

successful inflation compatible with the Planck constraint $r<0.058$ and
perturbative unitarity $\xi_{\sigma}<1$, which occurs for $5\times
10^{-13}\lesssim\tilde{\lambda}_{\sigma}\lesssim 5\times 10^{-10}$,

  * •

vacuum stability, which requires $|\lambda_{H\sigma}|\lesssim Y_{ii}^{2}$ and
$Y_{ii}\sim y\sim\lambda_{\sigma}^{1/4}$,

  * •

cold dark matter constituted by axions, which requires
$10^{10}\,\mathrm{GeV}\lesssim v_{\sigma}\lesssim 2.2\times
10^{11}\,\mathrm{GeV}$,111 The value of $f_{a}$ (or $v_{\sigma}$) accounting
for the observed cold dark matter abundance is subjected to the uncertainty in
the calculation of the relic axion abundance. Here we take the result of Ref.
[32] as a lower limit and that of Ref. [33] as an upper limit. See Sec. 5 for
more discussion.

the critical temperature of the PQ phase transition is predicted to be in the
range

$\displaystyle
T_{c}\sim\lambda_{\sigma}^{1/4}\,v_{\sigma}\sim\mathcal{O}(10^{7}\text{--}10^{9})\,\mathrm{GeV}.$
(2.19)

Around the critical temperature, there is a change in the equation of state of
the primordial plasma, which will be worked out precisely in the next section.

##  3 Equation of state for the Peccei-Quinn phase transition

In this section we review the computation of the effective numbers of
relativistic degrees of freedom, $g_{\ast\rho}(T)$ and $g_{\ast s}(T)$, which
affect the evolution of tensor perturbations after horizon crossing and thus
leave an imprint in the observable power spectrum, as in Eq. (1.5). The former
quantities are related to the equation of state of the thermal plasma. They
are linked to the total energy and entropy densities as

$\displaystyle\rho=\frac{\pi^{2}}{30}g_{\ast\rho}T^{4},\qquad
s=\frac{2\pi^{2}}{45}g_{\ast s}T^{3}.$ (3.1)

The first principle of thermodynamics and the Maxwell relation derived from
the enthalpy function allow to relate the former quantities to the system’s
pressure:

$\displaystyle\rho=T\frac{\partial p}{\partial T}-p,\qquad s=\frac{\partial
p}{\partial T}.$ (3.2)

In turn, for a system in equilibrium at constant temperature and volume and in
the absence of chemical potentials, the pressure can be related to minus the
Helmholtz free-energy per unit volume, which itself is related to the
partition function of the system at finite temperature. But the latter is also
connected to the effective potential:

$\displaystyle Z=e^{-\beta F}=e^{\beta Vp}=e^{-\beta VV_{\rm
eff}(T)}\Rightarrow p=-V_{\rm eff}(T),$ (3.3)

where $V$ is the three volume. As the system is assumed to be in equilibrium,
i.e. in the configuration that minimizes the free-energy, the potential
$V_{\rm eff}(T)$ is understood to be evaluated at its minimum. Thus, to
compute $g_{\ast\rho}$ and $g_{\ast s}$, one just needs to know the values of
the effective potential at its minimum as a function of temperature:

$\displaystyle\begin{aligned}
g_{\ast\rho}=&\,\frac{30}{\pi^{2}T^{4}}\left(V_{\rm
eff,min}(T)-T\frac{\partial V_{\rm eff,min}(T)}{\partial T}\right),\\\ g_{\ast
s}=&\,-\frac{45}{2\pi^{2}T^{3}}\frac{\partial V_{\rm eff,min}(T)}{\partial
T}.\end{aligned}$ (3.4)

In our calculations, we construct the effective potential in the Landau gauge
(with gauge-fixing parameter $\xi=0$) including the one-loop Coleman-Weinberg
potential, the one-loop thermal contributions, as well as higher-order QCD
corrections up to three-loop order. As we focus on high temperatures at which
the Higgs is stabilized at zero, we can consider a potential depending only on
$\sigma$:

$\displaystyle V_{\rm eff}(\sigma,T)=V(H=0,\sigma)+V^{\rm
CW}(\sigma,T)+V^{T}(\sigma,T)+V^{\rm QCD}(T).$ (3.5)

The tree-level piece $V(H=0,\sigma)$ follows from Eq. (2.1), while $V^{\rm
CW}(\sigma,T)$ designates the Coleman-Weinberg piece, $V^{T}(\sigma,T)$ is the
thermal potential, and $V^{\rm QCD}(T)$ accounts for QCD contributions to the
pressure. The temperature dependence in $V^{\rm CW}(\sigma,T)$ arises from the
use of an improved daisy resummation in which the effective masses for bosonic
fields appearing in $V^{\rm CW}$ (and $V^{T}$) are substituted by appropriate
temperature-dependent masses that incorporate the one-loop thermal corrections
to the corresponding self-energies at zero momentum. Such a resummation is
necessary because thermal corrections to the self-energies can dominate over
their zero-temperature counterparts, such that diagrams that differ by
insertions of thermal corrections to the propagators are all of the same order
and need to be resummed. To leading order in momenta, the corrections to the
propagators can be understood as thermal contributions to the masses. The
resummed diagrams are known as “daisy diagrams”, and in the usual daisy
resummation scheme, one simply incorporates the leading order thermal
corrections to the bosonic masses in a high-temperature expansion, which go as
$T^{2}$. However, as the high-temperature expansion is insensitive to mass
thresholds, the usual daisy resummation is not compatible with decoupling, so
that it cannot be used to follow the changes of $g_{\ast\rho}$ and $g_{\ast
s}$ across a phase transition, which are partly due to the decoupling of
particles from the plasma as the temperature drops below their mass. Hence it
becomes crucial to implement a resummation that is compatible with thermal
decoupling. In order to do so, we use an improved daisy resummation, which
resums the self-energies at zero momentum but keeps the full temperature
dependence.

For computing $V^{\rm CW}(\sigma,T)$ we use the Landau gauge result in the
$\overline{\rm MS}$ scheme,

$\displaystyle V^{\rm CW}(\sigma,T)=\frac{1}{64\pi^{2}}$
$\displaystyle\left[\sum_{V}m^{4}_{V}(\sigma,T)\left(\log\frac{m^{2}_{V}(\sigma,T)}{\mu^{2}}-\frac{5}{6}\right)+\sum_{S}m^{4}_{S}(\sigma,T)\left(\log\frac{m^{2}_{S}(\sigma,T)}{\mu^{2}}-\frac{3}{2}\right)\right.$
(3.6)
$\displaystyle\left.-\sum_{F}m^{4}_{F}(\sigma,T)\left(\log\frac{m^{2}_{F}(\sigma,T)}{\mu^{2}}-\frac{3}{2}\right)\right],$

where $V,S,F$ denote massive gauge bosons, real scalars, and Weyl fermions,
respectively, and where the sum over massive vectors includes a sum over the
three polarizations that propagate in the Landau gauge, while the sum over
fermion states also goes over the two spin/helicity states of each Weyl
fermion. $\mu$ designates the renormalization scale, while
$m^{2}_{V/S/F}(\sigma,T)$ are the field-dependent masses in the background of
$\sigma$, including the thermal corrections from the improved daisy
resummation detailed in Appendix A. Note that massless bosons and ghosts do
not contribute to the Coleman-Weinberg potential.

The one-loop thermal potential $V^{T}$ is given by

$\displaystyle
V^{T}=\frac{T^{4}}{2\pi^{2}}\left[\sum_{B}J_{B}\left(\frac{m^{2}_{B}(\sigma,T)}{T^{2}}\right)-\sum_{F}J_{F}\left(\frac{m^{2}_{F}(\sigma,T)}{T^{2}}\right)-\sum_{G}J_{B}\left(0\right)\right],$
(3.7)

where $B$ denotes all bosonic fields (massive and massless gauge bosons,
counting three polarizations, and real scalars), $F$ denotes all Weyl fermions
(counting 2 helicity states), and $G$ denotes complex ghosts—one for every
generator of every gauge group—which are massless in the Landau gauge. Note
that, in contrast to $V^{\rm CW}$, the massless fields contribute, which is
crucial to get the correct physical values of $g_{\ast\rho}$ and $g_{\ast s}$.
Despite ghosts being complex, there is no extra factor of 2 in the terms
corresponding to the ghosts because, even when the propagator in the Landau
gauge corresponds to 3 propagating polarizations, the fourth polarization
still contributes to the logarithm of the fluctuation determinant (after
appropriate subtraction of an infinite piece going as $\log\xi$), and the
corresponding term is cancelled by half of the ghost contribution. The
remaining ghost part and the three remaining gauge polarizations contribute as
in (3.7).222For example, given an Abelian field coupled to a complex scalar
one would sum over 3 gauge polarizations, 2 real scalars, and subtract one
ghost contribution; then effectively one has four degrees of freedom, as
expected in either a broken phase (three massive gauge polarizations and a
real scalar) or in the unbroken phase (2 massless gauge polarizations plus 2
real scalars). The functions $J_{B}$ and $J_{F}$ are given by

$\displaystyle\begin{aligned}
J_{B}(x)=\int_{0}^{\infty}dy\,y^{2}\log\left[1-\exp(-\sqrt{x+y^{2}})\right],\\\
J_{F}(x)=\int_{0}^{\infty}dy\,y^{2}\log\left[1+\exp(-\sqrt{x+y^{2}})\right].\end{aligned}$
(3.8)

In the ultrarelativistic limit, using the results
$J_{B}(0)=-(8/7)J_{F}(0)=-\pi^{4}/45$ in Eq. (3.7), one recovers the standard
pressure for a relativistic gas with the degrees of freedom in SMASH,

$\displaystyle-V_{T}\rightarrow p_{\rm
rel}=\frac{\pi^{2}}{90}\,g_{\ast\rho,\rm rel}T^{4},\quad g_{\ast\rho,\rm
rel}=124.5.$ (3.9)

The naive value of $g_{\ast\rho}$ above will be modified even in the high-
temperature limit as a consequence of the thermal corrections to the masses
and the higher-order QCD corrections. Regarding the daisy resummed masses, the
improved daisy resummation has to be applied to bosonic degrees of freedom
that couple to fields whose masses may become larger than the temperature. For
fields that remain light, the usual daisy resummation suffices. The PQ phase
transition gives masses to the real and imaginary components of $\sigma$, the
vector quark $Q$, and the right-handed neutrinos. The former couple to the
Higgs, the $\sigma$ fields, and the $SU(3)$ and hypercharge gauge bosons, for
which we should in principle use the improved daisy resummation. We will do so
for the Higgs, $\sigma$, and hypercharged fields, while for the $SU(3)$ gauge
bosons we will use an alternative treatment, to be discussed below, in order
to allow for the inclusion of higher-order QCD corrections without incurring
into double counting.

For a generic real scalar field $\phi_{j}$, the contributions of bosons and
fermions to its thermal mass beyond the high-temperature expansion for scalar
backgrounds $\phi_{j}=\bar{\phi}_{j}$ go as (see appendix A for details)

$\displaystyle\Delta
m^{2}_{\phi_{j}}(\bar{\phi}_{i},T)\supset\sum_{B}\,\frac{T^{2}}{\pi^{2}}J_{B}^{\prime}\left(\frac{m^{2}_{B}(\bar{\phi}_{i})}{T^{2}}\right)\left.\frac{\partial
m^{2}_{B}(\phi_{i})}{\partial\phi_{j}^{2}}\right|_{\phi_{i}\rightarrow\bar{\phi}_{i}}-2\sum_{F}\,\frac{T^{2}}{\pi^{2}}\,J_{F}^{\prime}\left(\frac{m^{2}_{F}(\bar{\phi}_{i})}{T^{2}}\right)\left.\frac{\partial
m^{2}_{F}(\phi_{i})}{\partial\phi_{j}^{2}}\right|_{\phi_{i}\rightarrow\bar{\phi}_{i}}.$
(3.10)

Again, the sum in $B$ goes over real scalars, and the sum in $F$ over Weyl
fermions. The masses $m^{2}_{B}(\phi_{i})$ and $m^{2}_{F}(\phi_{i})$ refer to
an arbitrary background for all the scalars at zero temperature, while
$m^{2}_{B}(\bar{\phi}_{i})$ and $m^{2}_{F}(\bar{\phi}_{i})$ refer to the
specific background of interest. In this way the factors $\partial
m^{2}_{B}(\phi_{i})/\partial{\phi_{j}}^{2}$ and $\partial
m^{2}_{F}(\phi_{i})/\partial{\phi_{j}}^{2}$ are just a convenient way to write
the portal couplings of $\phi_{j}$ to other scalars, or the Yukawa couplings
with fermions. Similarly, for an Abelian $U(1)$ field $G$ with gauge coupling
$\tilde{g}$, one has the following contribution to the effective mass due to
heavy Weyl fermions $F$ with charges $q_{F}$:

$\displaystyle\begin{aligned} \Delta
m^{2}_{G}(\bar{\phi}_{i},T)\supset\sum_{F}\frac{\tilde{g}^{2}q_{F}^{2}T^{2}}{\pi^{2}}K_{F}\left(\frac{m^{2}_{F}(\bar{\phi}_{i})}{T^{2}}\right),\end{aligned}$
(3.11)

where

$\displaystyle
K_{F}(x)=\int_{0}^{\infty}dy\,\frac{y^{2}e^{\sqrt{x+y^{2}}}}{(e^{\sqrt{x+y^{2}}}+1)^{2}}.$
(3.12)

For the contributions of fields that remain light in the range of temperatures
considered (e.g. the SM fields), one can use the standard daisy resummation
results.

As mentioned before one could also implement an improved daisy resummation for
the gluons, but this would complicate the inclusion of additional QCD loop
corrections without double counting. We have opted for including known three-
loop QCD corrections in the contribution $V^{\rm QCD}$, and implementing
decoupling of the heavy quarks $Q$ by interpolating between a seven flavour
regime for $T>m_{Q}$ (with $m_{Q}$ denoting the mass of the heavy quarks) and
a six flavour regime for $T<m_{Q}$. We use the results for $V^{\rm
QCD}(T)=-p^{\rm QCD}(T)$ up to order $g^{6}_{3}$ of Ref. [34], which account
for a variable number of massless flavours. The interpolation is performed by
considering a weighted sum of the seven and six flavour results. The weight of
the seven flavour contribution is taken as
$J_{F}(m^{2}_{Q}/T^{2})/J_{F}(0)$—which is zero at low temperature, and
approaches one at high $T$—while the weight of the six flavour contribution is
chosen as $1-J_{F}(m^{2}_{Q}/T^{2})/J_{F}(0)$. A consistent inclusion of the
three-loop QCD corrections requires knowing the QCD contributions to the beta
function of $g_{3}$ in SMASH up to three loops. For this we use the results of
Ref. [35], which give the following three-loop contribution to the beta
function $\beta_{g_{3}}=\mu\,\partial g_{3}/\partial\mu$ of SMASH:

$\displaystyle\Delta\beta_{g_{3}}=\frac{12629}{221184}\frac{g_{3}^{7}}{\pi^{6}}.$
(3.13)

With the full thermal potential $V_{\rm eff}(\sigma,T)$ of Eq. (3.7) computed
as explained above, one can then obtain the quantities $g_{\ast\rho}$ and
$g_{\ast s}$ from Eq. (3.4). However, a remaining subtlety is that the
formulae assume that all the relevant particles involved share the same
temperature. This neglects effects related to the loss of chemical
equilibrium. For massive particles that decouple from the thermal bath, even
though their effective temperature is expected to start to differ from that of
the rest of the plasma, the effect is not relevant for the computation of
$g_{\ast\rho}$ and $g_{\ast s}$, because decoupling ensures that the former
quantities are only sensitive to the light species. However, this is not so
for the axion degree of freedom emerging after the breaking of the PQ
symmetry. As corresponds to a pseudo-Goldstone, the axion interacts very
weakly, and at some point will fall out of chemical equilibrium with the rest
of the plasma. After that, the axion stops interacting with the thermal bath,
and the entropies of the former and the plasma are separately conserved,
implying a deviation of the temperature of the axion bath—which will be
referred to as $T_{\rm axion}$—from the temperature $T$ of the rest of the
plasma.

In order to estimate the corrections arising from the axion decoupling, one
needs to know the decoupling temperature $T_{\rm dec}$. In principle this can
be determined by computing thermally averaged rates or solving appropriate
Boltzmann equations,333Note that the setup considered here is somewhat
different from what usually assumed in the studies of the thermal axion
production such as Refs. [36, 37, 38]. Based on the results of these analyses
one may naively expect that axions decouple from the thermal plasma when the
temperature drops below the following critical value [37],  $\displaystyle
T_{\rm dec}\sim 1.7\times
10^{9}\,\mathrm{GeV}\,\left(\frac{f_{a}}{10^{11}\,\mathrm{GeV}}\right)^{2.246}.$
(3.14) However, in SMASH we cannot apply this estimate since the critical
temperature of the PQ phase transition is comparable or even lower than the
value of $T_{\rm dec}$ shown above. In other words, axions have not yet
emerged as a Goldstone degree of freedom at the would-be decoupling
temperature (3.14), and decouple only when the PQ symmetry breaking field
$\sigma$ acquires a sufficiently large expectation value. This fact motivates
us to consider $T_{\rm dec}$ to be lower than the estimate given by Eq (3.14),
as discussed in the text above Eq. (3.15). but here we adopt an approximation
consisting in estimating the temperature at which the phase transition can be
considered complete. While the usual critical temperature $T_{c}$ corresponds
to the point at which a new vacuum starts to arise, for a second-order phase
transition the vacuum at $T_{c}$ still lies at the origin, and the scale of
symmetry breaking only emerges at lower temperatures. An indicator of the
emergence of the PQ breaking scale is the trace of the stress-energy momentum,
which is zero in a plasma of massless particles but acquires peaks when new
dimensionful scales emerge and the associated massive degrees of freedom are
not yet thermally decoupled. Due to this, we will estimate $T_{\rm dec}$ as
the temperature at which the trace of the stress-energy tensor exhibits a
local maximum. For this we will consider the following dimensionless quantity
(called the trace anomaly),

$\displaystyle\Delta(T)=\frac{T^{\mu}_{\mu}}{T^{4}}=\frac{\rho-3p}{T^{4}}=\frac{1}{T^{4}}\left(4V_{\rm
eff,min}(T)-T\frac{\partial V_{\rm eff,min}(T)}{\partial T}\right).$ (3.15)

Once $T_{\rm dec}$ is determined from the maximum of $\Delta$, the values of
$g_{\ast\rho}$ and $g_{\ast s}$ accounting for the axion decoupling can be
calculated as follows. First, we may denote the quantities obtained from Eq.
(3.4) under the chemical equilibrium assumption as $g^{\rm eq}_{\ast\rho}$ and
$g^{\rm eq}_{\ast s}$. As a massless axion contributes one unit to both
$g^{\rm eq}_{\ast\rho}$ and $g^{\rm eq}_{\ast s}$, we can obtain the entropy
and energy density of the axion (with temperature $T_{\rm axion}$) and of the
thermal bath with the axion excluded (with temperature $T$) as:

$\displaystyle\begin{aligned} \rho^{\rm
axion}=&\,\frac{\pi^{2}}{30}\,T^{4}_{\rm axion},&\rho^{\rm
bath}=&\,\frac{\pi^{2}}{30}\,g_{\ast\rho}^{\rm bath}\,T^{4},&g^{\rm
bath}_{\ast\rho}=&\,g^{\rm eq}_{\ast\rho}-1,\\\ s^{\rm
axion}=&\,\frac{2\pi^{2}}{45}\,T^{3}_{\rm axion},&s^{\rm
bath}=&\,\frac{2\pi^{2}}{45}\,g_{\ast s}^{\rm bath}\,T^{3},&g^{\rm bath}_{\ast
s}=&\,g^{\rm eq}_{\ast s}-1.\end{aligned}$ (3.16)

Below the decoupling temperature $T_{\rm dec}$, the separate conservation of
$s^{\rm axion}$ and $s^{\rm bath}$ implies

$\displaystyle T_{\rm axion}=\left\\{\begin{array}[]{cl}T,&T\geq T_{\rm
dec},\\\ \left(\frac{g^{\rm bath}_{\ast s}(T)}{g^{\rm bath}_{\ast s}(T_{\rm
dec})}\right)^{\frac{1}{3}}T,&T<T_{\rm dec}.\end{array}\right.$ (3.19)

Then the final values of $g_{\ast\rho}$ and $g_{\ast s}$ follow from the total
energy density and entropy, $\rho=\rho^{\rm bath}+\rho^{\rm axion}$, $s=s^{\rm
bath}+s^{\rm axion}$, through Eqs. (3.1) and (3.16):

$\displaystyle g_{\ast\rho}=g^{\rm eq}_{\ast\rho}-1+\left(\frac{T_{\rm
axion}}{T}\right)^{4},\qquad g_{\ast s}=g^{\rm eq}_{\ast
s}-1+\left(\frac{T_{\rm axion}}{T}\right)^{3},$ (3.20)

where $T_{\rm axion}/T$ is to be computed using (3.19).

In Fig. 3, we plot the estimate of $g_{\ast\rho}$ and $g_{\ast s}$ obtained
based on the procedure described above. We also show the uncertainty of
$g_{\ast\rho}$ and $g_{\ast s}$ estimated by varying the renormalization scale
in the range $\mu\in(0.5\dots 2)m_{\rho}$, where
$m_{\rho}=\sqrt{2\lambda_{\sigma}}v_{\sigma}$ is the mass of the radial
direction $\rho$ of the PQ field, and the value of the unknown coefficient
$q_{c}$ in the order $g_{3}^{6}$ result for the QCD pressure [34] in the range
$q_{c}\in(-5000\dots+5000)$. The uncertainty from the $\mathcal{O}(g_{3}^{6})$
QCD pressure becomes less important at higher temperatures since the QCD
coupling becomes weaker and such a contribution is more suppressed at higher
energies. As shown in Fig. 3, the values of $g_{\ast\rho}$ and $g_{\ast s}$
start to fall off when the temperature becomes lower than the critical
temperature $T_{c}$ of the PQ phase transition.444The critical temperature of
the PQ phase transition shown in Figs. 3, 4 and Table 2 is defined as the
temperature at which the expectation value of the $\sigma$ field starts to
deviate from $\langle|\sigma|^{2}\rangle=0$, which is found numerically by
minimizing the full effective potential (3.5) rather than using the
approximate formula (2.18). Subsequently, there are step-like changes in
$g_{\ast\rho}$ and $g_{\ast s}$ when the temperature becomes about an order of
magnitude lower than $T_{c}$, which corresponds to the decoupling of heavy
fermions $Q$ and $N_{i}$. In the next section, we will see that such a change
in $g_{\ast\rho}$ and $g_{\ast s}$ leads to an observable signature in the
spectrum of primordial GWs.

Figure 3:  Temperature dependence of $g_{\ast\rho}(T)$ (top left) and $g_{\ast
s}(T)$ (top right) obtained from Eq. (3.20). Shaded regions represent the
uncertainty due to the choice of the renormalization scale and estimate of the
QCD corrections (see text). Note that these regions do not necessary indicate
the largest source of uncertainty, since we expect that there might be further
corrections if we go beyond the sudden decoupling approximation to estimate
the contribution of relativistic axions. Gray dashed lines correspond to the
critical temperature of the PQ phase transition. In the bottom panel, central
values of $g_{\ast\rho}$ and $g_{\ast s}$ are compared directly. In these
figures, parameter values are taken according to benchmark point 1 in Table 3.

Practically, the calculation based on the thermal potential (3.5) exhibits
numerical instabilities at low temperatures, since it involves cancellation of
large numbers. This is clearly seen in Fig. 3, where the values of
$g_{\ast\rho}$ and $g_{\ast s}$ fluctuate a lot at
$T=\mathcal{O}(10^{5})\,\mathrm{GeV}$. This fact prevents us from tracing the
decoupling of the radial component of the PQ field with mass $m_{\rho}\sim
10^{5}$–$10^{7}\,\mathrm{GeV}$ and evaluating the effective relativistic
degrees of freedom at lower temperatures. To avoid this issue, we evaluate
thermodynamic quantities at lower temperatures by adding the contribution of
the free $\rho$ particle, which can be estimated analytically, to $\Delta$ for
the SM:

$\displaystyle\Delta(T)$ $\displaystyle=\Delta_{\rm SM}(T)+\Delta_{\rho}(T),$
(3.21) $\displaystyle\Delta_{\rho}(T)$
$\displaystyle=\left.\frac{x}{\pi^{2}}J_{B}^{\prime}(x)\right|_{x\,=\,m_{\rho}^{2}/T^{2}},$
(3.22)

where $\Delta_{\rm SM}$ denotes the contribution of the SM particles. Here we
use the result of Ref. [5] to evaluate $\Delta_{\rm SM}$ and match the SMASH
result (3.15) to Eq. (3.21). As shown in Fig. 4, the two results indeed agree
with each other, which allows us to interpolate between the SM and SMASH
results by using Eq. (3.21). The interpolation is performed by applying a
Gaussian filter to smooth out the fluctuations in the SMASH result (3.15) at
low temperatures and switching from Eq. (3.15) to Eq. (3.21) at a point where
the central values of two results cross each other. After obtaining the
interpolated function of $\Delta(T)$, we estimate the pressure $p(T)$ by
integrating $\Delta(T)$ in a similar way to Ref. [5], to obtain $g_{\ast\rho}$
and $g_{\ast s}$ at arbitrary temperature.

Figure 4:  Temperature dependence of $\Delta(T)$ evaluated according to Eq.
(3.21) (blue) and Eq. (3.15) (orange). Right panel just shows a zoomed plot of
the left panel. Shaded regions represent the uncertainty due to the choice of
the renormalization scale and estimate of the QCD corrections (see text), and
the gray dashed line corresponds to the critical temperature of the PQ phase
transition. In these figures, parameter values are taken according to
benchmark point 1 in Table 3.

To finalize, we provide four SMASH benchmark points associated with different
values of the effective inflationary quartic coupling at high scales. These
imply different effective quartic couplings $\tilde{\lambda}_{\sigma}$ along
the inflationary valley (see Eq. (2.5)). Under the assumption of a given value
of the non-minimal gravitational coupling $\xi_{\sigma}$, and under the
requirement of a consistent reheating history, the value of the tensor-to-
scalar ratio $r(k_{*})$ during inflation becomes fixed in terms of
$\tilde{\lambda}_{\sigma}$ [12, 13]. The chosen benchmark points are
summarized in Table 3. In all four cases, we take the scale of the PQ symmetry
breaking to be $v_{\sigma}=1.2\times 10^{11}\,\mathrm{GeV}$. In Table 2, we
also show the values of the critical temperature of the PQ phase transition
$T_{c}$ and the axion decoupling temperature $T_{\rm dec}$, which corresponds
to the local maximum of $\Delta$, for each benchmark point.

Table 1: Parameter values for the chosen SMASH benchmark points. We take $v_{\sigma}=1.2\times 10^{11}\,\mathrm{GeV}$ in all four benchmark points. $\phi_{*}$ denotes the value of the inflaton field (in the Jordan frame) when the Planck pivot scale of 0.002 ${\rm Mpc}^{-1}$ crosses the horizon. Benchmark point | $T_{c}\,[10^{9}\,\mathrm{GeV}]$ | $T_{\rm dec}\,[10^{7}\,\mathrm{GeV}]$
---|---|---
1 | 0.48–0.51 | 3.9–4.1
2 | 1.25–1.31 | 8.7–9.1
3 | 1.34–1.44 | 5.6–5.8
4 | 4.50–5.20 | 15.6–16.1
Table 2: Values of the critical temperature of the PQ phase transition $T_{c}$
and axion decoupling temperature $T_{\rm dec}$ for the chosen SMASH benchmark
points. The range of values corresponds to the uncertainty due to the choice
of the renormalization scale.

In Fig. 3, we show the estimate of $g_{\ast\rho}$, $T_{\rm axion}/T$, and
$\Delta$ for the chosen SMASH benchmark points. We see that the value of the
effective relativistic degrees of freedom in SMASH becomes larger than that in
the SM at high temperatures. In particular, there are two step-like features
at $T\sim m_{\rho}$ (threshold of $\rho$ particle) and at $T\sim T_{\rm dec}$
(threshold of heavy fermions $Q$ and $N_{i}$). These events correspond to the
two peaks in $\Delta(T)$, which are shown in the bottom panel of Fig. 3. Note
that the ordering of the benchmark points (and hence the ordering of the value
of the tensor-to-scalar ratio $r$) does not imply an analogous ordering of the
locations of the peaks, since $r$ is related to the effective inflationary
coupling $\tilde{\lambda}_{\sigma}$, which itself depends on a combination of
three coupling parameters (see Eq. (2.5)). Furthermore, the overall values of
$g_{\ast\rho}$ and $g_{\ast s}$ remain slightly larger than the SM values even
at lower temperatures, because of the extra contribution from relativistic
axions (see Eq. (3.20)). The contribution of the relativistic axions is
proportional to a power of $T_{\rm axion}/T$, which decays at low temperature
like $\propto[g_{\ast s}^{\rm bath}(T)]^{1/3}$ as shown in the top right panel
of Fig. 3.

Figure 5: Temperature dependence of $g_{\ast\rho}(T)$ (top left), $T_{\rm
axion}/T$ (top right), and $\Delta(T)$ (bottom) for the chosen SMASH benchmark
points. For the plots of $g_{\ast\rho}(T)$, we also show the SM result (gray
dashed line) and its uncertainty (gray shaded region) estimated in Ref. [5].

Before closing this section, let us estimate the asymptotic values of
$g_{\ast\rho}$ and $g_{\ast s}$ after the neutrino decoupling, which are
relevant to the present observations. The conservation of entropy implies that
the ratio of the effective temperatures of neutrinos and axions to photon
temperature after the neutrino decoupling are respectively

$\displaystyle\frac{T_{\nu}}{T_{\gamma}}=\left(\frac{4}{11}\right)^{\frac{1}{3}},\quad\frac{T_{\rm
axion}}{T_{\gamma}}=\left(\frac{43}{11g_{\ast s}^{\rm bath}(T_{\rm
dec})}\right)^{\frac{1}{3}}.$ (3.39)

By using these formulae, the asymptotic values of $g_{\ast\rho}$ and $g_{\ast
s}$ after the neutrino decoupling can be estimated as

$\displaystyle g_{\ast\rho,0}$
$\displaystyle=2+\frac{21}{4}\left(\frac{4}{11}\right)^{\frac{4}{3}}+\left(\frac{43}{11g_{\ast
s}^{\rm bath}(T_{\rm dec})}\right)^{\frac{4}{3}}\simeq
3.36+0.01\,\left(\frac{114}{g_{\ast s}^{\rm bath}(T_{\rm
dec})}\right)^{\frac{4}{3}},$ (3.40) $\displaystyle g_{\ast s,0}$
$\displaystyle=\frac{43}{11}\left(1+\frac{1}{g_{\ast s}^{\rm bath}(T_{\rm
dec})}\right)\simeq 3.91+0.03\,\left(\frac{114}{g_{\ast s}^{\rm bath}(T_{\rm
dec})}\right).$ (3.41)

The former can be recast in terms of an excess in effective number of
neutrinos,

$\displaystyle\Delta N_{\rm eff}=\frac{4}{7}\left(\frac{T_{\rm
axion}}{T_{\nu}}\right)^{4}\simeq 0.0245\,\left(\frac{114}{g_{\ast s}^{\rm
bath}(T_{\rm dec})}\right)^{\frac{4}{3}}.$ (3.42)

The excess in $N_{\rm eff}$ may be probed in future CMB and large scale
structure observations [39]. We numerically find that the chosen SMASH
benchmark points predict $g_{\ast s}^{\rm bath}(T_{\rm dec})\simeq 114$, and
hence $g_{\ast\rho,0}\simeq 3.37$, $g_{\ast s,0}\simeq 3.94$, and $\Delta
N_{\rm eff}\simeq 0.0245$. These values are slightly larger than those
predicted in the SM, $g_{\ast\rho,0}\simeq 3.36$, $g_{\ast s,0}\simeq 3.91$,
and $\Delta N_{\rm eff}=0$.555Strictly speaking, when we take account of the
feeble interaction of neutrinos with the electromagnetic plasma and the
leading order quantum electrodynamics corrections, the asymptotic values of
the effective relativistic degrees of freedom in the SM become
$g_{\ast\rho,0}\simeq 3.38$ and $g_{\ast s,0}\simeq 3.93$ [5] (see also Refs.
[40, 41, 42] for a precision calculation of the neutrino decoupling and Refs.
[43, 44, 45, 46] for more recent works), which are slightly larger than the
commonly used values $g_{\ast\rho,0}\simeq 3.36$ and $g_{\ast s,0}\simeq
3.91$. The SMASH predictions shown in Eqs. (3.40)–(3.42) may also be subjected
to similar corrections, whose precise estimation is beyond the scope of this
paper. As the amplitude of GWs observed today depends on $g_{\ast s,0}$ (see
Eq. (1.5)), we use the value estimated based on Eq. (3.41) when we evaluate
the GW signatures in the next section.

##  4 Spectrum of gravitational waves

Now we compute the spectrum of GWs predicted in SMASH by using the ingredients
obtained in the previous sections. The GWs can be described by spatial metric
perturbations $h_{ij}$ on the Friedmann-Robertson-Walker (FRW) background,

$\displaystyle ds^{2}=-dt^{2}+a^{2}(t)(\delta_{ij}+h_{ij})dx^{i}dx^{j},$ (4.1)

where it is assumed that $h_{ij}$ satisfy the transverse-traceless conditions,
$h^{i}_{i}=\partial^{i}h_{ij}=0$. We expand $h_{ij}$ in terms of their Fourier
modes,

$\displaystyle h_{ij}(t,{\bf
x})=\sum_{\lambda}\int\frac{d^{3}k}{(2\pi)^{3}}h^{\lambda}(t,{\bf
k})\epsilon^{\lambda}_{ij}({\bf k})e^{i{\bf k\cdot x}},$ (4.2)

where $\lambda=+,\times$ represents two independent polarization states, and
$\epsilon^{\lambda}_{ij}({\bf k})$ are the spin-2 polarization tensors that
satisfy
$\sum_{ij}\epsilon^{\lambda}_{ij}(\epsilon^{\lambda^{\prime}}_{ij})^{*}=2\delta^{\lambda\lambda^{\prime}}$.
In the analysis of the spectrum of the primordial GWs, it will be convenient
to introduce a dimensionless quantity $X(t,k)$ defined as

$\displaystyle h^{\lambda}(t,{\bf k})=h^{\lambda}_{{\bf k},p}X(t,k),$ (4.3)

where $h^{\lambda}_{{\bf k},p}$ represents the amplitudes of primordial tensor
perturbations that are fixed when the corresponding modes cross outside the
horizon, and $X(t,k)$ satisfies $X(t,k)\to 1$ for $k\ll aH$. The former is
related to the primordial tensor power spectrum in Eq. (1.1) as

$\displaystyle\langle h_{p,ij}(t,{\bf x})h_{p,ij}(t,{\bf
x})\rangle=\int\frac{dk}{k}\Delta_{t}^{2}(k),$ (4.4)

where $h_{p,ij}(t,{\bf x})$ are the primordial tensor fields given by Eq.
(4.2) with $h^{\lambda}(t,{\bf k})$ replaced by $h^{\lambda}_{{\bf k},p}$, and
$\langle\dots\rangle$ denotes the ensemble average. It can also be written as

$\displaystyle\Delta_{t}^{2}(k)=\frac{2k^{3}|h_{p}(k)|^{2}}{\pi^{2}},$ (4.5)

where $|h_{p}(k)|^{2}$ is defined by

$\displaystyle\langle h^{\lambda*}_{{\bf k},p}h^{\lambda^{\prime}}_{{\bf
k^{\prime}},p}\rangle=(2\pi)^{3}\delta^{\lambda\lambda^{\prime}}\delta^{(3)}({\bf
k}-{\bf k^{\prime}})|h_{p}(k)|^{2}.$ (4.6)

On the other hand, the quantity $X(t,k)$ describes the evolution of GWs after
inflation, and it can be well-approximated by the WKB solution $X\propto
a^{-1}e^{\pm ik\tau}$ for $k\gg aH$.

Substituting Eq. (4.2) into the formula for the energy density of GWs [47],

$\displaystyle\rho_{\rm gw}(t)=\frac{1}{32\pi G}\langle\dot{h}_{ij}(t,{\bf
x})\dot{h}_{ij}(t,{\bf x})\rangle,$ (4.7)

where $G$ is the Newton’s gravitational constant, we see that the spectrum of
GWs $\Omega_{\rm gw}(f)$ can be represented as the product of the primordial
tensor power spectrum $\Delta_{t}^{2}(f)$ and the transfer function
$\mathcal{T}_{0}(f)$ as shown in Eq. (1.2). Here, the transfer function is
calculated from a derivative of the dimensionless quantity $X(\tau,k)$ with
respect to conformal time $d\tau=dt/a$,

$\displaystyle\mathcal{T}_{0}(f)=\frac{1}{12a_{0}^{2}H_{0}^{2}}\left[\frac{dX(\tau_{0},k)}{d\tau}\right]^{2},$
(4.8)

where the comoving wavenumber $k$ is related to the frequency via $f=k/(2\pi
a_{0})$.

At linear order in perturbation theory, the evolution of $X(\tau,k)$ is
described by the following integro-differential equation [48, 49],

$\displaystyle\frac{d^{2}X(u)}{du^{2}}+\frac{2}{a(u)}\frac{da(u)}{du}\frac{dX(u)}{du}+X(u)=-24\sum_{i=\gamma,\nu,a}F_{i}(u)\left[\frac{1}{a(u)}\frac{da(u)}{du}\right]^{2}\int^{u}_{u_{i}}dU\left[\frac{j_{2}(u-U)}{(u-U)^{2}}\right]\frac{dX(U)}{dU},$
(4.9)

where $u=k\tau$ and $j_{n}(z)$ is the spherical Bessel function of the first
kind. The right-hand side of Eq. (4.9) represents the collisionless damping
effect due to free-streaming particles [48], which is proportional to the
fraction of the energy density in these particles,

$\displaystyle F_{i}(u)\equiv\frac{\rho_{i}(u)}{\rho_{\rm crit}(u)}.$ (4.10)

Here we consider photons ($i=\gamma$), neutrinos ($i=\nu$), and relativistic
axions ($i=a$) as the sources of the damping effect. The lower limit of the
integration $u_{i}$ is taken as the time at which the corresponding particles
decouple from the thermal bath. Following Ref. [5], we take the times
corresponding to $T=3000\,K$ and $T=2\,\mathrm{MeV}$ for photons and
neutrinos, respectively. For axions, we use $T=T_{\rm dec}$, which is
determined from the local maximum of $\Delta$ defined in Eq. (3.15). It turned
out that the contribution of free-streaming neutrinos leads to the suppression
of the amplitude of GWs by $\lesssim 35\,\%$ at frequencies
$10^{-16}\,\mathrm{Hz}\lesssim f\lesssim 10^{-10}\,\mathrm{Hz}$, and that of
photons leads to the additional damping by $\lesssim 14\,\%$ at $f\sim
10^{-17}\,\mathrm{Hz}$ [5]. In SMASH, there is an extra contribution from
relativistic axions, which we study in this section.666The effect of
relativistic axions on the spectrum of primordial GWs in a general context was
discussed in Refs. [50, 51].

We numerically solve Eq. (4.9) with initial conditions $X(0)=1$ and
$dX(0)/du=0$, together with the Friedmann equation that describes the
evolution of the scale factor $a(u)$. The numerical solution is evaluated up
to a finite value of $u=u_{\rm end}=20$ and extrapolated until the present
time $\tau_{0}$ by matching it to the WKB solution $X\propto a^{-1}e^{\pm
iu}$, which is substituted to Eq. (4.8). We confirmed that the value of
$\mathcal{T}_{0}h^{2}$ can be estimated with an accuracy of $\lesssim 0.2\,\%$
in our numerical scheme. After obtaining the transfer function, the spectrum
of GWs can be estimated by multiplying it by the primordial tensor power
spectrum as shown in Eq. (1.2).

In Fig. 6, we show the spectrum of GWs for parameters corresponding to the
benchmark point 1 in Table 3. As shown in this figure, various cosmological
events occurring at temperature $T=T_{\rm hc}$ are imprinted on the shape of
the GW spectrum at the corresponding frequency (cf. Eq. (1.4)). In this plot,
the feature due to the PQ phase transition, which we are interested in, is
less clear compared to that from $e^{+}e^{-}$ annihilation and the QCD
crossover, since the change in the relativistic degrees of freedom is milder
than in the latter cases. However, we will see below that such a feature could
be distinguishable if a GW detector has enough sensitivity.

Figure 6: The spectrum of GWs predicted in SMASH for a broad frequency
interval. The values of model parameters are fixed according to the benchmark
point 1 in Table 3. We also show the temperature $T_{\rm hc}$ corresponding to
the horizon crossing of the mode with frequency $f$ evaluated according to Eq.
(1.4).

The plot shown in Fig. 6 takes account of the fact that the amplitude of GWs
coherently oscillate with a phase $2k\tau_{0}$. Practically, such oscillations
cannot be resolved in direct detection experiments that are sensitive to GWs
at high frequencies, and we may estimate the amplitude of GWs by replacing a
rapidly oscillating factor by $1/2$. In the following, we use the notation
$\Omega_{\rm gw}$ (or $\mathcal{T}_{0}$) to denote the averaged quantity.

As described in Eq. (4.9), we have included the damping effect due to free-
streaming axions in our analysis of the GW spectrum. The coefficient
$F_{a}(u)$ defined by Eq. (4.10) represents the efficiency of the damping
effect, and it can be rewritten as

$\displaystyle F_{a}(u)=\frac{\frac{1}{2}\left(\frac{T_{\rm
axion}}{T}\right)^{4}\left(\frac{g_{\ast s,0}}{g_{\ast
s}(T)}\right)^{\frac{4}{3}}\Omega_{\gamma}h^{2}}{\Omega_{M}h^{2}\left(\frac{a(u)}{a_{0}}\right)+\frac{g_{\ast\rho}(T)}{2}\left(\frac{g_{\ast
s,0}}{g_{\ast s}(T)}\right)^{\frac{4}{3}}\Omega_{\gamma}h^{2}},$ (4.11)

where $\Omega_{M}h^{2}\simeq 0.14$ and $\Omega_{\gamma}h^{2}\simeq 2.47\times
10^{-5}$ are the present fraction of the energy densities of matter and
photons, respectively. In the radiation-dominated era, this function can be
simplified as

$\displaystyle F_{a}(u)\to\frac{\left(T_{\rm
axion}/T\right)^{4}}{g_{\ast\rho}(T)}.$ (4.12)

Figure 7 shows the impact of the axion damping effect on the spectrum of GWs.
Note that $F_{a}$ decreases monotonically for $T<T_{\rm dec}$, since
$F_{a}=(T_{\rm axion}/T)^{4}(g_{\ast\rho})^{-1}\propto(g^{\rm bath}_{\ast
s})^{4/3}(g_{\ast\rho})^{-1}\propto(g^{\rm bath}_{\ast s})^{1/3}$ for $g^{\rm
bath}_{\ast s}\simeq g_{\ast\rho}$. On the other hand, $F_{a}$ loses its
physical meaning for $T>T_{\rm dec}$ since axions do not exist as Goldstone
bosons in this regime. Therefore, the fraction of the energy density in
relativistic axions is at most only $\sim 0.8\,\%$, which leads to a $\lesssim
1\,\%$ suppression of the amplitude of GWs at frequencies $f\lesssim
1\,\mathrm{Hz}$. Although the effect is tiny, it might be relevant to future
high-sensitivity GW experiments, as we discuss below.

Figure 7: Effect of free-streaming axions on the spectrum of GWs. In the left
panel, the ratio of $\Omega_{\rm gw}h^{2}$ obtained by including the axion
damping effect to that obtained without including it (i.e. $F_{a}(u)=0$) is
plotted. In the right panel, $F_{a}(u)$ given by Eq. (4.11) is plotted as a
function of temperature $T$ or frequency $f$ that corresponds to $T_{\rm
hc}=T$. The shaded region in the right panel corresponds to $T>T_{\rm dec}$.
These plots are produced based on the model parameters corresponding to the
benchmark point 1 in Table 3.

Figure 8 shows the transfer function of GWs for the modes reentering the
horizon around the epoch of the PQ phase transition. It clearly shows that
there is a step-like feature at $f\sim 1\,\mathrm{Hz}$ corresponding to the
change in the effective relativistic degrees of freedom shown in Fig. 3. In
addition to this feature, the amplitude of GWs below that frequency is
slightly suppressed due to the damping effect from free-streaming axions.
There is also a minor step-like feature at $f\sim 10^{-2}\,\mathrm{Hz}$, which
is induced by the decoupling of $\rho$ particles. These features are
contrasted with the prediction in the SM, where the transfer function remains
almost flat at these frequency ranges.777The transfer function in the SM is
not exactly flat and slightly red-tilted, because of the renormalization group
running of gauge and Yukawa couplings [5]. The reason why the SMASH result
with no axion damping (orange dashed line in Fig. 8) is larger than the SM
(gray dotted line) at $f\lesssim 10^{-3}\,\mathrm{Hz}$ is just because the
transfer function is proportional to $g_{\ast s,0}^{4/3}$ (see Eq. (1.5)) and
its value in SMASH is larger than that in the SM, as shown in Eq. (3.41).

Figure 8:  Frequency dependence of the transfer function
$\mathcal{T}_{0}(f)h^{2}$ predicted in SMASH (blue solid line) and the SM
(gray dotted line). The result obtained by using the equation of state in
SMASH, but not including the damping effect from free-streaming axions (i.e.
$F_{a}(u)=0$), is also plotted as the orange dashed line for the sake of
comparison. The SMASH results shown in this figure are produced based on the
model parameters corresponding to the benchmark point 1 in Table 3.

The spectrum of GWs predicted in SMASH is compared with the projected
sensitivities of future GW experiments in Fig. 9. The experimental
sensitivities in this figure represent the power-law integrated sensitivity
(PLIS) curves [52] with a detection threshold given by the signal-to-noise
ratio (SNR) of $\mathrm{SNR}=1$, which implies that any power-law signal that
comes above these curves results in $\mathrm{SNR}>1$. The ongoing ground-based
experiments such as the Advanced Laser Interferometer Gravitational Wave
Observatory (aLIGO) [53], Advanced Virgo (aVirgo) [54], and the Kamioka
Gravitational Wave Detector (KAGRA) [55] may not be sensitive enough to probe
the range of GW amplitudes shown in this figure, while more advanced detectors
such as the Cosmic Explorer (CE) [56, 57] and the Einstein Telescope (ET) [58]
will improve the sensitivities at
$f\sim\mathcal{O}(1\text{--}10)\,\mathrm{Hz}$. For lower frequencies, the
planned space-borne experiment LISA [59, 60] is expected to probe GWs in the
range $f\gtrsim 1\,\mathrm{mHz}$, and future pulser timing array experiments,
the International Pulsar Timing Array (IPTA) [61] and the Square Kilometre
Array (SKA) [62, 63] will be sensitive to GWs at $f\sim 10^{-9}\,\mathrm{Hz}$.
Unfortunately, the nearly scale-invariant GW spectra predicted in the standard
slow-roll inflationary models like SMASH are beyond the reach of the
experiments mentioned above, but there is a possibility that they can be
probed by more advanced space-borne detectors like BBO [6, 7, 8, 9] and DECIGO
[10, 11]. These experiments are designed to probe GWs at the frequency range
$f\gtrsim 0.1\,\mathrm{Hz}$, which avoids the expected astrophysical
foregrounds coming from extragalactic WD binaries [64].888 In order to make a
fair comparison between the PLIS curves and WD confusion noise, we have used
the following equation to draw the WD curve in Fig. 9,
$\displaystyle\Omega_{\rm
WD}(f)=\frac{4\pi^{2}}{3H_{0}^{2}}\frac{f^{3}}{\sqrt{2\Delta t_{\rm obs}\Delta
f}}\frac{S_{\rm WD}(f)}{\Gamma(f)},$ (4.13) where $\Delta t_{\rm obs}$ is the
observing time, $\Delta f$ is a range of frequencies at which the integrand
giving the SNR becomes approximately constant, $S_{\rm WD}(f)$ is the strain
noise power spectrum for the WD foregrounds, and $\Gamma(f)$ is the overlap
reduction function of detectors. Equation (4.13) may represent a typical
amplitude of the WD confusion noise in the correlation analysis between two GW
detectors. Here we take $\Delta t_{\rm obs}=1\,\mathrm{year}$, $\Delta
f=f/10$, and $\Gamma(f)=3/5$, assuming observations in LISA/DECIGO-like
detectors, and adopt the fitting formula for $S_{\rm WD}(f)$ taken from Ref.
[65]. It is notable that the GW signals predicted by some of the benchmark
points in SMASH are accessible to DECIGO and BBO. On the other hand, in order
to probe the detailed features on the GW spectrum, further improvement in the
experimental sensitivities is indispensable. This may be achievable in an
idealized experiment whose sensitivity is limited by the quantum noise
(ultimate DECIGO) [10]. Intriguingly, the frequency range $f\sim
1\,\mathrm{Hz}$ corresponding to the best sensitivity range of ultimate DECIGO
almost exactly coincides with the range where the non-trivial feature due to
the PQ phase transition in SMASH emerges.

Figure 9:  Spectrum of GWs predicted by the chosen SMASH benchmark points
shown in Table 3 (thick solid lines) and projected sensitivities of future GW
experiments represented as PLIS curves with $\mathrm{SNR}=1$ (dashed lines).
All the sensitivity curves other than ultimate DECIGO are taken from Ref. [66]
(see also Ref. [67]). The sensitivity curve for ultimate DECIGO is estimated
by using the parameters specified in Ref. [68]. An estimate of the
astrophysical foregrounds from extragalactic WD binaries is also shown as gray
thin sold line (see footnote 8 for details).

More detailed comparisons between the GW signal in SMASH and the sensitivity
of ultimate DECIGO are shown in Fig. 4. In this figure, we have estimated the
uncertainty of $\Omega_{\rm gw}$ in ultimate DECIGO in terms of a discrete set
of frequency bins, following the procedure described in Appendix D of Ref.
[5].999The estimation of the sensitivity of ultimate DECIGO in this paper is
the same as Ref. [5] except for the following two modifications: First, we
take the angular efficiency factor (or the overall normalization of the
overlap reduction function) as $F=3/5$, which corresponds to the assumption
that the detectors consist of interferometers with opening angle of $\pi/3$
and make use of two independent data streams [67], rather than $F=2/5$, which
was used in Ref. [5] but corresponds to the case of interferometers with
perpendicular arms and a single data channel. Second, instead of assuming the
constant overlap reduction function (i.e. $\Gamma(f)=\text{const.}$) at all
frequencies, we include the decay of $\Gamma(f)$ at higher frequencies
according to Fig. 10 of Ref. [67]. The first modification leads to an overall
improvement of the sensitivity in $\Omega_{\rm gw}$ by a factor of $2/3$,
while the second one reduces the sensitivity at $f\gtrsim 10\,\mathrm{Hz}$. If
the amplitude of the primordial tensor perturbation is sufficiently large,
which is the case for the benchmark point 1 (top left panel of Fig. 4), the
primordial GW spectrum is more tilted because of the slow-roll consistency
relation (see Eq. (2.11)), and the feature due to the change in the effective
relativistic degrees of freedom in SMASH is seen as an additional suppression
of the tilted spectrum. In this case, the shape of the GW spectrum can be
identified over a wide frequency range, since the amplitude of GW itself is
large compared to the sensitivity of ultimate DECIGO. On the other hand, the
step-like feature becomes more manifest in the cases with smaller values of
$r$, but in such cases the experimental noise is more severe because of the
reduction of the overall GW amplitude. The worst case is shown in the bottom
right panel of Fig. 4 that corresponds to the benchmark point 4. In spite of
this difficulty, it is clear that in all four cases the GW signatures
resulting from the transfer function in SMASH are distinguishable from those
in the SM, which opens up the possibility of probing the non-trivial thermal
history predicted in SMASH at future high-sensitivity GW experiments.

Figure 10: The spectrum of primordial GWs predicted in SMASH benchmark point 1
(top left), 2 (top right), 3 (bottom left), and 4 (bottom right),
corresponding to the choices of parameter values shown in Table 3. Gray dotted
lines represent the results obtained by using the effective relativistic
degrees of freedom in the SM. Dark yellow bars correspond to the sensitivity
of ultimate DECIGO.

##  5 Discussion and conclusions

In this paper we have focused on the primordial GWs predicted in SMASH. We
have worked out precisely the unique prediction for the shape of the spectrum
of GWs due to the nontrivial feature of the PQ phase transition described in
this model. The thermodynamic quantities in SMASH across the second order PQ
phase transition were computed for four different benchmark parameter values
shown in Table 3 by adopting an improved daisy resummation scheme compatible
with thermal decoupling, which allowed us to estimate the effective
relativistic degrees of freedom $g_{\ast\rho}(T)$ and $g_{\ast s}(T)$ for a
wide temperature interval, as shown in Fig. 3. Using these estimates of
$g_{\ast\rho}(T)$ and $g_{\ast s}(T)$, we computed numerically the spectrum of
the primordial GWs originating from tensor fluctuations during inflation and
compared it to the sensitivity of ultimate DECIGO. It turned out that the
change in the relativistic degrees of freedom caused by the PQ phase
transition in SMASH results in a step-like feature in the spectrum of
primordial GWs, and that ultimate DECIGO is capable to identify such a signal
as shown in Fig. 4.

As shown in Eq. (1.4), the modifications in the spectrum of the primordial GWs
at particular frequency ranges can be related to cosmological events occurring
at corresponding temperatures. In this sense, the GW spectrum represents not
only a snapshot of the early universe but also a movie looking deep into the
very early universe, cf. Fig. 6. Different GW detectors sensitive to different
frequency ranges can give snapshots of different epochs, providing us with a
possibility to reconstruct this “cosmic movie” by combining the results of
various GW observations. The SMASH model can be seen as a benchmark model for
such a strategy, as it is a minimal model of particle cosmology, allowing us
to describe the whole thermal history of the universe from inflation until
today in a decisive way as explicitly worked out in this paper. Furthermore,
one can also think of other minimal models of particle cosmology, such as for
example the Neutrino Minimal SM ($\nu$MSM) [69, 70], a model based on a $B-L$
symmetry [71]101010Attempts to estimate the spectrum of GWs predicted in the
model proposed in Ref. [71] were made in Refs. [72, 73]., and a model based on
a flavour-dependent global $U(1)$ symmetry [74], and study the spectrum of GWs
predicted in such models. The formalism developed in this paper may also be
used for such endeavors.

Given that real direct observations of the primordial GWs in far advanced
detectors like ultimate DECIGO will become feasible within this century [10],
an important and necessary step in the near future is that the next generation
CMB experiments should find primordial B-mode polarization patterns and
measure an exact value of the tensor-to-scalar ratio $r$ at the CMB scale. If
it will be the case, we will be able to resolve the uncertainty regarding the
value of $r$ (and hence the value of the effective inflationary coupling
$\tilde{\lambda}_{\sigma}$) and make more decisive predictions for the
spectrum of GWs in SMASH at scales relevant to the space-borne GW experiments.
Then, at the time when such direct detection experiments start operating, one
can obtain a richer information about the shape of the spectrum of primordial
GWs and hence about viable parameter values in SMASH, by combining the results
of CMB and direct detection experiments. Such joint measurements of the GW
spectrum over broad frequency ranges in various different detectors were
considered recently in Ref. [75].

Throughout this paper, we have assumed a certain value for the scale of the PQ
symmetry breaking or axion decay constant, $v_{\sigma}=f_{a}=1.2\times
10^{11}\,\mathrm{GeV}$. In axion models with the post-inflationary PQ symmetry
breaking and domain wall number $N_{\rm DW}=1$ (including SMASH), the value of
$f_{a}$ should be determined from the requirement that the axion becomes the
main constituent of dark matter. Unfortunately, the determination of $f_{a}$
based on this cosmological argument currently involves a large uncertainty due
to the fact that the calculation of the axion dark matter abundance from
global strings remains far from straightforward, despite a lot of effort made
by several groups [76, 77, 78, 79, 33, 80, 81, 82, 83, 84] (see the third
bullet point between Eqs. (2.18) and (2.19)). The latest results of numerical
simulations of axionic strings reported in Ref. [32], which take account of
not only the radiation of axions in the scaling regime but also the dynamics
at the nonlinear regime when the effect of the axion mass becomes relevant,
imply that $f_{a}\lesssim 10^{10}\,\mathrm{GeV}$ when one extrapolates the
simulation results to the realistic value of the string tension. On the other
hand, the analysis in Ref. [33] based on the effective action method that
directly simulates strings with large tension [79] gives a larger value,
$f_{a}=(2.21\pm 0.29)\times 10^{11}\,\mathrm{GeV}$. A careful assessment to
understand the source of discrepancy between these results is still missing at
the moment.

The determination of the value of $f_{a}$ is crucial for the calculation of
the spectrum of the primordial GWs in SMASH, since it affects the critical
temperature of the PQ phase transition (see Eqs. (2.18) and (2.19)) and
changes the frequency range at which the nontrivial feature on the GW spectrum
emerges. We generically expect that the signal from the PQ phase transition
may become more difficult to detect if the required value of $f_{a}$ becomes
lower, since the feature arising from the change in the effective relativistic
degrees of freedom will shift to lower frequencies and be hidden by the WD
foreground contamination noise appearing at $f\lesssim 0.1\,\mathrm{Hz}$.

Independently of the theoretical predictions of $f_{a}$ mentioned above, its
value will be fixed definitely once the axion dark matter is found in the
future haloscope experiments, such as the Axion Dark Matter eXperiment (ADMX)
[85], Capp Ultra Low Temperature Axion Search in Korea (CULTASK) [86],
Haloscope at Yale Sensitive to Axion CDM Experiment (HAYSTAC) [87], MAgnetized
Disc and Mirror Axion eXperiment (MADMAX) [88, 89], and Relic Axion Detector
Exploratory Setup (RADES) [90]. It is likely that the axion dark matter would
be detected already at the epoch when the advanced GW detectors like ultimate
DECIGO are taking data, if it actually constitutes the main component of dark
matter. Hence, the measurements of the spectrum of GWs should bring us some
further information about the model, such as the self coupling of the PQ
field, the Higgs portal coupling, and the number and masses of exotic
fermions, in addition to the value of the PQ scale. In this sense, the future
GW observations can be used to probe the details of the PQ sector that may not
be reached in other high-energy experiments.

## Acknowledgments

We would like to thank Valerie Domcke for discussions and useful comments. AR
acknowledges partial funding by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) under Germany’s Excellence Strategy – EXC 2121
“Quantum Universe”. Part of this research was supported by the Munich
Institute for Astro- and Particle Physics (MIAPP) which is funded by DFG under
Germany’s Excellence Strategy – EXC-2094 – 390783311. This work of KS is
supported by Leading Initiative for Excellent Young Researchers (LEADER), the
Ministry of Education, Culture, Sports, Science, and Technology (MEXT), Japan.
CT acknowledges financial support by SFB 1258 of DFG.

##  Appendix A Improved daisy resummation compatible with thermal decoupling

As the values of $g_{\ast\rho}$ and $g_{\ast s}$ depend on temperature, an
accurate calculation has to incorporate all relevant thermal effects. As is
well-known from thermal field theory, thermal corrections for two-point
functions can become of the same order as their zero-temperature counterparts.
Indeed, in the high-temperature limit and at low momentum it can be argued
that thermal contributions to bosonic self-energies are of the order of
$g_{i}^{2}T^{2}/m^{2}_{\rm tree}$ times the zero-temperature result, where
$g_{i}$ stands for couplings of the theory, and $m_{\rm tree}$ is a zero-
temperature mass. Therefore, a resummation is needed for $g_{i}T\gtrsim m_{\rm
tree}$. In particular, this becomes necessary for temperatures above and near
the critical temperature of phase transitions, as the latter are associated
with effective masses of the form

$\displaystyle m^{2}_{X,\rm eff}(T)=m^{2}_{X,\rm tree}+\Sigma_{X}(p=0,T)\sim
m^{2}_{X,\rm tree}+f_{X}(g_{i})T^{2}(1+O(m^{2}_{X,\rm tree}/T^{2}))\,,$ (A.1)

becoming zero for some particular bosonic field $X$ with self-energy
$\Sigma_{X}$. The resummation is achieved by substituting two-point functions
by their thermally corrected counterparts, which is equivalent to resumming
one-loop self-energy insertions in so-called “daisy” diagrams, as illustrated
in Fig. 11. In the low momentum limit, the resummation of self-energy
insertions can be implemented as a thermal shift of the masses of the
different bosonic particles as

$\displaystyle m^{2}_{X,\rm tree}\rightarrow m^{2}_{X,\rm
tree}+\Sigma_{X}(p=0,T)=m^{2}_{X,\rm eff}(T).$ (A.2)

Figure 11: The upper two lines illustrate the corrections to the bosonic
propagators coming from the scalars and fermions that become heavy during the
PQ phase transition. The Daisy resummation in the effective potential accounts
for all one-loop diagrams with external background legs and arbitrary
insertions of the corrections to the internal two-point functions. This is
illustrated at the bottom of the figure, where the dots in the external lines
represent the background fields.

Traditionally, one simply employs the leading corrections to
$\Sigma_{X}(p=0,T)$ in a high-temperature expansion, as in Eq. (A.1). However,
this approach is doomed to fail as the temperature goes down and the high-
temperature expansion loses accuracy. This becomes crucially problematic for
estimating the variation of the effective numbers of relativistic degrees of
freedom, $g_{\ast\rho}$ and $g_{\ast s}$, across a phase transition. The
former quantities change during the transition precisely because some
particles in the thermal plasma become nonrelativistic, such that $T/m_{X,\rm
eff}$ drops below one and the particles decouple from the plasma. This
decoupling is reflected by Boltzmann suppression factors in the corresponding
contributions to the self-energies, which are not captured with the
traditional daisy resummation. Hence, if one uses it to estimate
$g_{\ast\rho}$ and $g_{\ast s}$ the results become unphysical.

Nevertheless, one can improve the usual daisy resummation by avoiding the
high-$T$ expansion in the thermal corrections to the self-energies, while
still keeping the zero-momentum approximation. As will be seen next, the
resulting thermal corrections to the self-energies (which are absorbed into
corrections to the masses) can be then expressed in terms of thermal integrals
related to the finite-temperature corrections to the effective potential (e.g.
the functions (3.8) appearing in (3.7)). We first consider the contribution to
the thermal mass of a scalar singlet induced by other scalars and fermions, to
be followed by the calculation of the contribution to the mass of an Abelian
gauge boson induced by fermions. This covers the relevant cases in SMASH,
where the fields becoming heavy during the phase transition are the heavy
quarks and the real component of the complex scalar singlet. These particles
contribute to the thermal masses of the Higgs boson, the hypercharge gauge
boson, and the complex scalar of SMASH.

For the first case we assume a real scalar $\phi$ with a Yukawa interaction
with a Dirac fermion $\psi$ and a quartic interaction with a real scalar
$\chi$:

$\displaystyle{\cal
L}\supset\left(-m_{\psi}\bar{\psi}\psi-\frac{\hat{y}}{\sqrt{2}}\phi\overline{\psi}P_{L}\psi+c.c.\right)-\frac{m^{2}_{\phi}}{2}\phi^{2}-\frac{m^{2}_{\chi}}{2}\chi^{2}-\frac{\lambda_{\phi}}{4!}\,\phi^{4}-\frac{\lambda_{\chi}}{4!}\,\chi^{4}-\frac{\lambda}{4}\,\phi^{2}\chi^{2}.$
(A.3)

For applications to SMASH, one can either interpret $\chi$ as a real component
of the Higgs doublet and $\phi$ as a real component of the complex singlet, or
the other way around; one simply has to use different mappings of the
couplings in SMASH to the couplings in Eq. (A.3). In the following, we will
assume for simplicity that the Yukawa coupling $\hat{y}$ is real. For future
reference, the thermal contributions to the effective potential of $\phi$ (we
assume zero background for $\chi$) can be separated as:

$\displaystyle\begin{aligned} V^{T}=&\,V^{\phi,T}+V^{\chi,T}+V^{\psi,T},\\\
V^{\phi/\chi,T}=&\,\frac{T^{4}}{2\pi^{2}}\,J_{B}\left(\frac{m^{2}_{\phi/\chi}(\phi)}{T^{2}}\right),\\\
V^{\psi,T}=&\,-4\frac{T^{4}}{2\pi^{2}}\,J_{F}\left(\frac{m^{2}_{\psi}(\phi)}{T^{2}}\right).\end{aligned}$
(A.4)

Note the extra factor of 2 in $V^{\psi,T}$, as we are dealing with a Dirac
fermion composed of 2 Weyl fermions. The effective masses in a $\phi$
background are

$\displaystyle\begin{aligned}
m^{2}_{\phi}(\phi)=&\,m^{2}_{\phi}+\frac{\lambda_{\phi}}{2}\phi^{2},\\\
m^{2}_{\chi}(\phi)=&\,m^{2}_{\chi}+\frac{\lambda}{2}\phi^{2},\\\
m_{\psi}(\phi)=&\,m_{\psi}+\frac{\hat{y}}{\sqrt{2}}\phi.\end{aligned}$ (A.5)

We have the following one-loop corrections at finite $T$ for the $\phi$ self-
energy in a $\phi$ background:

$\displaystyle\begin{aligned}
-i\Sigma_{\phi}(0)=&\,-i\Sigma^{\phi}_{\phi}(0)-i\Sigma^{\chi}_{\phi}(0)-i\Sigma^{\psi}_{\phi}(0),\\\
-i\Sigma^{\phi}_{\phi}(0)=&\,\frac{i\lambda_{\phi}}{2\beta}\sum_{n=-\infty}^{\infty}\int\frac{d^{3}l}{(2\pi)^{3}}\frac{1}{(2\pi
n)^{2}-\vec{l}^{2}-m^{2}_{\chi}(\phi)},\\\
-i\Sigma^{\chi}_{\phi}(0)=&\,\frac{i\lambda}{2\beta}\sum_{n=-\infty}^{\infty}\int\frac{d^{3}l}{(2\pi)^{3}}\frac{1}{(2\pi
n)^{2}-\vec{l}^{2}-m^{2}_{\chi}(\phi)},\\\
-i\Sigma^{\psi}_{\phi}(0)=&\,-\frac{i\hat{y}^{2}}{2\beta}\sum_{n=-\infty}^{\infty}\int\frac{d^{3}l}{(2\pi)^{3}}{\rm
Tr}\,\frac{\not{l}\not{l}+O(m^{2}_{\psi}(\phi))}{\left((2\pi(n+1/2))^{2}-\vec{l}^{2}-m^{2}_{\psi}(\phi)\right)^{2}},\\\
=&\,-\frac{2i\hat{y}^{2}}{\beta}\sum_{n=-\infty}^{\infty}\int\frac{d^{3}l}{(2\pi)^{3}}\frac{1}{(2\pi(n+1/2))^{2}-\vec{l}^{2}-m^{2}_{\psi}(\phi)}+O(\hat{y}^{2}m^{2}_{\psi}(\phi)),\end{aligned}$
(A.6)

where we neglected contributions proportional $\hat{y}^{2}m^{2}_{\phi}(\phi)$.
The reason is that we focus on realizations with no tree-level mass terms for
the fermions, so that $m^{2}_{\phi}(\phi)\propto\hat{y}^{2}$, and thus the
terms proportional to $\hat{y}^{2}m^{2}_{\phi}(\phi)$ are $O(\hat{y}^{4})$ and
can be safely neglected for $\hat{y}\ll 1$. Using the identity (see e.g. [91])

$\displaystyle\begin{aligned}
\frac{1}{\beta}\sum_{n=-\infty}^{\infty}f(p^{0}=i\omega_{n})=&\,\int_{-i\infty}^{i\infty}\frac{dz}{4\pi
i}(f(z)+f(-z))+\eta_{F/B}\int_{-i\infty+\epsilon}^{i\infty+\epsilon}\frac{dz}{2\pi
i}n_{F/B}(z)(f(z)+f(-z)),\\\ \eta_{F}=&\,-1,\quad\eta_{B}=1,\quad
n_{F/B}(z)=\frac{1}{e^{\beta z}-\eta_{F/B}},\end{aligned}$ (A.7)

one can separate the loop integrals in Eq. (A.6) into zero $T$ and finite $T$
contributions, corresponding respectively to the first and second terms in Eq.
(A.7). Starting with the $\Sigma^{\chi}_{\phi}$ contribution, keeping only the
finite-$T$ terms leads to

$\displaystyle\begin{aligned}
-i\Sigma^{\chi,T}_{\phi}(0)=&\,\frac{\lambda}{2\pi}\int_{-i\infty+\epsilon}^{i\infty+\epsilon}dz\int\frac{d^{3}l}{(2\pi)^{3}}\frac{1}{e^{\beta
z}-1}\cdot\frac{1}{z^{2}-\vec{l}^{2}-m^{2}_{\chi}(\phi)}.\end{aligned}$ (A.8)

Closing the integration contour of $z$ on the ${\rm Re}\,z>0$ hemiplane, and
picking up the residue at $z=E_{l}=\sqrt{\vec{l}^{2}+m^{2}_{\chi}(\phi)}$,
with

$\displaystyle{\text{Res}}\,\left.\frac{1}{e^{\beta
z}-1}\cdot\frac{1}{z^{2}+e^{2}_{l}}\right|_{z=E_{l}}=\frac{1}{2E_{l}}\frac{1}{e^{\beta
E_{l}}-1},$ (A.9)

we get

$\displaystyle\begin{aligned}
\Sigma^{\chi,T}_{\phi}(0)=&\,\frac{\lambda}{2}\int\frac{d^{3}l}{(2\pi)^{3}}\frac{1}{E_{l}}\frac{1}{e^{\beta
E_{l}}-1}=\frac{\lambda}{4\pi^{2}}\int
l^{2}dl\frac{1}{\sqrt{l^{2}+m^{2}_{\chi}(\phi)}}\frac{1}{e^{\beta\sqrt{l^{2}+m^{2}_{\chi}(\phi)}}-1}.\end{aligned}$
(A.10)

Changing variables to $y=l/T$, and defining $x\equiv
m^{2}_{\chi}(\phi)/T^{2}$, we arrive at

$\displaystyle\begin{aligned}
\Sigma^{\chi,T}_{\phi}(0)&\,=\frac{T^{2}\lambda}{4\pi^{2}}\int
dy\frac{y^{2}}{\sqrt{x+y^{2}}}\frac{1}{e^{\beta\sqrt{x+y^{2}}}-1}.\end{aligned}$
(A.11)

From Eq. (3.8) we see that this is related to derivatives of the thermal loop
functions, and in fact to derivatives of the corresponding thermal
contribution to the effective potential

$\displaystyle\begin{aligned}
\Sigma^{\chi,T}_{\phi}(0)=&\,\frac{T^{2}\lambda}{2\pi^{2}}J_{B}^{\prime}\left(\frac{m^{2}_{\chi}(\phi)}{T^{2}}\right)=\frac{T^{2}}{\pi^{2}}J_{B}^{\prime}\left(\frac{m^{2}_{\chi}(\phi)}{T^{2}}\right)\frac{dm^{2}_{\chi}(\phi)}{d\phi^{2}}=2\frac{d}{d\phi^{2}}V^{\chi,T}.\end{aligned}$
(A.12)

Of course, the last identity just confirms the fact that
$\Sigma^{\chi,T}_{\phi}(0)$ can be interpreted as a contribution to the mass
of $\phi$, since the total mass goes $m^{2}\sim
d^{2}V/(d^{2}\phi)|_{\phi=0}=2dV/(d\phi^{2})|_{\phi=0}$.

With the analytic expansion of $J_{B}$,

$\displaystyle J_{B}(x)=-\frac{\pi^{4}}{45}+\frac{\pi^{2}}{12}x+O(x^{3/2}),$
(A.13)

we recover the shift of the mass of $\phi$ due to thermal $\chi$ loops in the
high temperature limit,

$\displaystyle\Delta^{\chi}m^{2}_{\phi}=\Sigma_{\phi}^{\chi,T}(0)=\frac{\lambda}{24}\,T^{2}+\dots\,,$
(A.14)

which is the result employed in traditional daisy resummations. As noted
before, this is incompatible with decoupling, while the result that avoids the
high $T$ expansion is general and perfectly compatible with decoupling:

$\displaystyle\Delta^{\chi}m^{2}_{\phi}=\frac{T^{2}}{\pi^{2}}J_{B}^{\prime}\left(\frac{m^{2}_{\chi}(\phi)}{T^{2}}\right)\frac{dm^{2}_{\chi}(\phi)}{d\phi^{2}}.$
(A.15)

The same reasoning applies to the contribution $\Sigma^{\phi,T}_{\phi}$. In
the case of the fermions, repeating the same arguments as above, one has

$\displaystyle\Delta^{\psi}m^{2}_{\phi}=-\frac{4T^{2}}{\pi^{2}}J_{F}^{\prime}\left(\frac{m^{2}_{\psi}(\phi)}{T^{2}}\right)\frac{dm^{2}_{\psi}(\phi)}{d\phi^{2}}=2\frac{d}{d\phi^{2}}V^{\psi,T}.$
(A.16)

With the analytic expansion of $J_{F}$,

$\displaystyle J_{F}(x)=\frac{7\pi^{4}}{360}-\frac{\pi^{2}}{24}x+O(x^{3/2}),$
(A.17)

one recovers the usual result for daisy resummations

$\displaystyle\Delta^{\psi}m^{2}_{\phi}=\frac{\hat{y}^{2}}{12}\,T^{2}+\dots.$
(A.18)

Finally, we can consider the thermal mass of an Abelian boson $G$ (with
associated gauge coupling $g$) induced by a heavy Dirac fermion $\psi$ with
charge $q$. In the low momentum limit the thermal mass is captured by minus
the temporal component of the self-energy [92],

$\displaystyle\Delta^{\psi}m^{2}_{G}=-\Pi_{G}^{\psi,00}(0)=-\frac{q^{2}g^{2}}{\beta}\sum_{n=-\infty}^{\infty}\int\frac{d^{3}l}{(2\pi)^{3}}\,{\rm
Tr}\,\frac{(\not{l}+m_{\psi})\gamma^{0}(\not{l}+m_{\psi})\gamma^{0}}{\left((2\pi(n+1/2))^{2}-\vec{l}^{2}-m^{2}_{\psi}\right)^{2}}.$
(A.19)

Proceeding in analogous manner as before the result is

$\displaystyle\Delta^{\psi}m^{2}_{G}=\frac{2g^{2}q^{2}T^{2}}{\pi^{2}}K_{F}\left(\frac{m^{2}_{\psi}}{T^{2}}\right),$
(A.20)

with $K_{F}$ given by Eq. (3.12). In the relativistic limit one has
$K_{F}(0)={\pi^{2}/6}$, giving the standard daisy resummation

$\displaystyle\Delta^{\psi}m^{2}_{G}=\frac{g^{2}q^{2}}{3}\,T^{2}+\dots.$
(A.21)

### Summary

The previous formulae can be easily generalized to arbitrary backgrounds
involving more than one real scalar component, and for higher numbers of
scalars and fermions. Labelling real scalars as $\phi_{i}$, and considering a
background $\phi_{i}=\bar{\phi}_{i}$, the improved daisy resummation of self-
energy corrections of a real scalar $\phi_{j}$ can be implemented by
considering the following shift in the mass of $\phi_{j}$, arising from heavy
real scalars $B=\\{\phi_{i}\\}$ and Weyl fermions $F$:

$\displaystyle\Delta
m^{2}_{\phi_{j}}(\bar{\phi}_{i},T)\supset\sum_{B}\,\frac{T^{2}}{\pi^{2}}J_{B}^{\prime}\left(\frac{m^{2}_{B}(\bar{\phi}_{i})}{T^{2}}\right)\left.\frac{\partial
m^{2}_{B}(\phi_{i})}{\partial\phi_{j}^{2}}\right|_{\phi_{i}\rightarrow\bar{\phi}_{i}}-2\sum_{F}\,\frac{T^{2}}{\pi^{2}}\,J_{F}^{\prime}\left(\frac{m^{2}_{F}(\bar{\phi}_{i})}{T^{2}}\right)\left.\frac{\partial
m^{2}_{F}(\phi_{i})}{\partial\phi_{j}^{2}}\right|_{\phi_{i}\rightarrow\bar{\phi}_{i}}.$
(A.22)

The derivatives of the thermal loop functions implement decoupling, as they go
to zero exponentially as $m^{2}_{X}/T\gg 1$. One can use numerical
interpolating functions in the code. For the fields coupling to $\phi$ and
remaining light (e.g. additional gauge bosons, ignored in the above formula),
one can use the shift in the mass arising in the standard daisy resummation.

For the improved daisy resummation of the self-energy of an Abelian gauge
field coupling to heavy Weyl fermions $F$, one can instead use the following
shift in the mass:

$\displaystyle\Delta
m^{2}_{G}\supset\sum_{F}\frac{g^{2}q^{2}T^{2}}{\pi^{2}}K_{F}\left(\frac{m^{2}_{G}(\bar{\phi}_{i})}{T^{2}}\right).$
(A.23)

The function $K_{F}$ becomes again exponentially suppressed if the argument is
large, as expected from decoupling. For additional light fields coupling to
the Abelian boson one can follow the standard treatment.

## References

  * [1] L. P. Grishchuk, “Amplification of gravitational waves in an isotropic universe,” Sov. Phys. JETP 40 (1975) 409 [Zh. Eksp. Teor. Fiz. 67 (1974) 825].
  * [2] A. A. Starobinsky, “Spectrum of relict gravitational radiation and the early state of the universe,” JETP Lett. 30 (1979) 682 [Pisma Zh. Eksp. Teor. Fiz. 30 (1979) 719].
  * [3] D. J. Schwarz, “Evolution of gravitational waves through cosmological transitions,” Mod. Phys. Lett. A 13 (1998) 2771 [gr-qc/9709027].
  * [4] N. Seto and J. Yokoyama, “Probing the equation of state of the early universe with a space laser interferometer,” J. Phys. Soc. Jap. 72 (2003), 3082-3086 [arXiv:gr-qc/0305096 [gr-qc]].
  * [5] K. Saikawa and S. Shirai, “Primordial gravitational waves, precisely: The role of thermodynamics in the Standard Model,” JCAP 1805 (2018) 035 [arXiv:1803.01038 [hep-ph]].
  * [6] S. Phinney et al., The Big Bang Observer: Direct detection of gravitational waves from the birth of the Universe to the Present, NASA Mission Concept Study (2004).
  * [7] J. Crowder and N. J. Cornish, “Beyond LISA: Exploring future gravitational wave missions,” Phys. Rev. D 72 (2005) 083005 [gr-qc/0506015].
  * [8] V. Corbin and N. J. Cornish, “Detecting the cosmic gravitational wave background with the big bang observer,” Class. Quant. Grav. 23 (2006) 2435 [gr-qc/0512039].
  * [9] G. M. Harry, P. Fritschel, D. A. Shaddock, W. Folkner and E. S. Phinney, “Laser interferometry for the big bang observer,” Class. Quant. Grav. 23 (2006) 4887 Erratum: [Class. Quant. Grav. 23 (2006) 7361].
  * [10] N. Seto, S. Kawamura and T. Nakamura, “Possibility of direct measurement of the acceleration of the universe using 0.1-Hz band laser interferometer gravitational wave antenna in space,” Phys. Rev. Lett. 87 (2001) 221103 [astro-ph/0108011].
  * [11] S. Kawamura et al., “The Japanese space gravitational wave antenna DECIGO,” Class. Quant. Grav. 23 (2006) S125.
  * [12] G. Ballesteros, J. Redondo, A. Ringwald and C. Tamarit, “Unifying inflation with the axion, dark matter, baryogenesis and the seesaw mechanism,” Phys. Rev. Lett. 118 (2017) no.7, 071802 [arXiv:1608.05414 [hep-ph]].
  * [13] G. Ballesteros, J. Redondo, A. Ringwald and C. Tamarit, “Standard Model—axion—seesaw—Higgs portal inflation. Five problems of particle physics and cosmology solved in one stroke,” JCAP 1708 (2017) 001 [arXiv:1610.01639 [hep-ph]].
  * [14] G. Ballesteros, J. Redondo, A. Ringwald and C. Tamarit, “Several Problems in Particle Physics and Cosmology Solved in One SMASH,” Front. Astron. Space Sci. 6 (2019) 55 [arXiv:1904.05594 [hep-ph]].
  * [15] K. N. Abazajian et al. [CMB-S4 Collaboration], “CMB-S4 Science Book, First Edition,” arXiv:1610.02743 [astro-ph.CO].
  * [16] T. Matsumura et al., “Mission design of LiteBIRD,” J. Low. Temp. Phys. 176 (2014) 733 [arXiv:1311.2847 [astro-ph.IM]].
  * [17] P. Ade et al. [Simons Observatory Collaboration], “The Simons Observatory: Science goals and forecasts,” JCAP 1902 (2019) 056 [arXiv:1808.07445 [astro-ph.CO]].
  * [18] R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38 (1977) 1440.
  * [19] S. Weinberg, “A New Light Boson?,” Phys. Rev. Lett. 40 (1978) 223.
  * [20] F. Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40 (1978) 279.
  * [21] J. Preskill, M. B. Wise and F. Wilczek, “Cosmology of the Invisible Axion,” Phys. Lett. B 120 (1983) 127.
  * [22] L. F. Abbott and P. Sikivie, “A Cosmological Bound on the Invisible Axion,” Phys. Lett. B 120 (1983) 133.
  * [23] M. Dine and W. Fischler, “The Not So Harmless Axion,” Phys. Lett. B 120 (1983) 137.
  * [24] P. Minkowski, “$\mu\to e\gamma$ at a Rate of One Out of $10^{9}$ Muon Decays?,” Phys. Lett. B 67 (1977) 421.
  * [25] M. Gell-Mann, P. Ramond and R. Slansky, “Complex Spinors and Unified Theories,” Conf. Proc. C 790927 (1979) 315 [arXiv:1306.4669 [hep-th]].
  * [26] T. Yanagida, “Horizontal Symmetry And Masses Of Neutrinos,” Conf. Proc. C 7902131 (1979) 95 [Conf. Proc. C 7902131 (1979) 95].
  * [27] R. N. Mohapatra and G. Senjanovic, “Neutrino Mass and Spontaneous Parity Violation,” Phys. Rev. Lett. 44 (1980) 912.
  * [28] M. Fukugita and T. Yanagida, “Baryogenesis Without Grand Unification,” Phys. Lett. B 174 (1986) 45.
  * [29] G. Ballesteros and J. A. Casas, “Large tensor-to-scalar ratio and running of the scalar spectral index with Instep Inflation,” Phys. Rev. D 91 (2015) 043502 [arXiv:1406.3342 [astro-ph.CO]].
  * [30] N. Aghanim et al. [Planck], “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020), A6 [arXiv:1807.06209 [astro-ph.CO]].
  * [31] Y. Akrami et al. [Planck], “Planck 2018 results. X. Constraints on inflation,” Astron. Astrophys. 641 (2020), A10 [arXiv:1807.06211 [astro-ph.CO]].
  * [32] M. Gorghetto, E. Hardy and G. Villadoro, “More Axions from Strings,” [arXiv:2007.04990 [hep-ph]].
  * [33] V. B. Klaer and G. D. Moore, “The dark-matter axion mass,” JCAP 11 (2017), 049 [arXiv:1708.07521 [hep-ph]].
  * [34] K. Kajantie, M. Laine, K. Rummukainen and Y. Schroder, “The Pressure of hot QCD up to $g^{6}\ln(1/g)$,” Phys. Rev. D 67 (2003) 105008 [hep-ph/0211321].
  * [35] A. G. M. Pickering, J. A. Gracey and D. R. T. Jones, “Three loop gauge beta function for the most general single gauge coupling theory,” Phys. Lett. B 510 (2001), 347-354 [arXiv:hep-ph/0104247 [hep-ph]].
  * [36] E. Masso, F. Rota and G. Zsembinszki, “On axion thermalization in the early universe,” Phys. Rev. D 66 (2002), 023004 [arXiv:hep-ph/0203221 [hep-ph]].
  * [37] P. Graf and F. D. Steffen, “Thermal axion production in the primordial quark-gluon plasma,” Phys. Rev. D 83 (2011), 075011 [arXiv:1008.4528 [hep-ph]].
  * [38] A. Salvio, A. Strumia and W. Xue, “Thermal axion production,” JCAP 01 (2014), 011 [arXiv:1310.6982 [hep-ph]].
  * [39] D. Baumann, D. Green and B. Wallisch, “Searching for light relics with large-scale structure,” JCAP 08 (2018), 029 [arXiv:1712.08067 [astro-ph.CO]].
  * [40] G. Mangano, G. Miele, S. Pastor and M. Peloso, “A Precision calculation of the effective number of cosmological neutrinos,” Phys. Lett. B 534 (2002), 8-16 [arXiv:astro-ph/0111408 [astro-ph]].
  * [41] G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti and P. D. Serpico, “Relic neutrino decoupling including flavor oscillations,” Nucl. Phys. B 729 (2005), 221-234 [arXiv:hep-ph/0506164 [hep-ph]].
  * [42] P. F. de Salas and S. Pastor, “Relic neutrino decoupling with flavour oscillations revisited,” JCAP 07 (2016), 051 [arXiv:1606.06986 [hep-ph]].
  * [43] J. J. Bennett, G. Buldgen, M. Drewes and Y. Y. Y. Wong, “Towards a precision calculation of the effective number of neutrinos $N_{\rm eff}$ in the Standard Model I: The QED equation of state,” JCAP 03 (2020), 003 [arXiv:1911.04504 [hep-ph]].
  * [44] M. Escudero Abenza, “Precision early universe thermodynamics made simple: $N_{\rm eff}$ and neutrino decoupling in the Standard Model and beyond,” JCAP 05 (2020), 048 [arXiv:2001.04466 [hep-ph]].
  * [45] K. Akita and M. Yamaguchi, “A precision calculation of relic neutrino decoupling,” JCAP 08 (2020), 012 [arXiv:2005.07047 [hep-ph]].
  * [46] J. Froustey, C. Pitrou and M. C. Volpe, “Neutrino decoupling including flavour oscillations and primordial nucleosynthesis,” [arXiv:2008.01074 [hep-ph]].
  * [47] M. Maggiore, “Gravitational Waves. Vol. 1: Theory and Experiments,” Oxford University Press, Oxford U.K. (2007).
  * [48] S. Weinberg, “Damping of tensor modes in cosmology,” Phys. Rev. D 69 (2004), 023503 [arXiv:astro-ph/0306304 [astro-ph]].
  * [49] Y. Watanabe and E. Komatsu, “Improved Calculation of the Primordial Gravitational Wave Spectrum in the Standard Model,” Phys. Rev. D 73 (2006), 123515 [arXiv:astro-ph/0604176 [astro-ph]].
  * [50] R. Jinno, T. Moroi and K. Nakayama, “Probing dark radiation with inflationary gravitational waves,” Phys. Rev. D 86 (2012), 123502 [arXiv:1208.0184 [astro-ph.CO]].
  * [51] R. Jinno, T. Moroi and K. Nakayama, “Inflationary Gravitational Waves and the Evolution of the Early Universe,” JCAP 01 (2014), 040 [arXiv:1307.3010 [hep-ph]].
  * [52] E. Thrane and J. D. Romano, “Sensitivity curves for searches for gravitational-wave backgrounds,” Phys. Rev. D 88 (2013) no.12, 124032 [arXiv:1310.5300 [astro-ph.IM]].
  * [53] J. Aasi et al. [LIGO Scientific], “Advanced LIGO,” Class. Quant. Grav. 32 (2015), 074001 [arXiv:1411.4547 [gr-qc]].
  * [54] F. Acernese et al. [VIRGO], “Advanced Virgo: a second-generation interferometric gravitational wave detector,” Class. Quant. Grav. 32 (2015) no.2, 024001 [arXiv:1408.3978 [gr-qc]].
  * [55] T. Akutsu et al. [KAGRA], “KAGRA: 2.5 Generation Interferometric Gravitational Wave Detector,” Nature Astron. 3 (2019) no.1, 35-40 [arXiv:1811.08079 [gr-qc]].
  * [56] B. P. Abbott et al. [LIGO Scientific], “Exploring the Sensitivity of Next Generation Gravitational Wave Detectors,” Class. Quant. Grav. 34 (2017) no.4, 044001 [arXiv:1607.08697 [astro-ph.IM]].
  * [57] D. Reitze et al., “Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO,” Bull. Am. Astron. Soc. 51, 035 [arXiv:1907.04833 [astro-ph.IM]].
  * [58] M. Punturo et al., “The Einstein Telescope: A third-generation gravitational wave observatory,” Class. Quant. Grav. 27 (2010), 194002.
  * [59] P. Amaro-Seoane et al. [LISA], “Laser Interferometer Space Antenna,” [arXiv:1702.00786 [astro-ph.IM]].
  * [60] J. Baker et al., “The Laser Interferometer Space Antenna: Unveiling the Millihertz Gravitational Wave Sky,” [arXiv:1907.06482 [astro-ph.IM]].
  * [61] G. Hobbs et al., “The international pulsar timing array project: using pulsars as a gravitational wave detector,” Class. Quant. Grav. 27 (2010), 084013 [arXiv:0911.5206 [astro-ph.SR]].
  * [62] G. Janssen et al., “Gravitational wave astronomy with the SKA,” PoS AASKA14 (2015), 037 [arXiv:1501.00127 [astro-ph.IM]].
  * [63] A. Weltman et al., “Fundamental physics with the Square Kilometre Array,” Publ. Astron. Soc. Austral. 37 (2020), e002 [arXiv:1810.02680 [astro-ph.CO]].
  * [64] A. J. Farmer and E. S. Phinney, “The gravitational wave background from cosmological compact binaries,” Mon. Not. Roy. Astron. Soc. 346 (2003), 1197 [arXiv:astro-ph/0304393 [astro-ph]].
  * [65] A. Nishizawa, K. Yagi, A. Taruya and T. Tanaka, “Cosmology with space-based gravitational-wave detectors — dark energy and primordial gravitational waves —,” Phys. Rev. D 85 (2012), 044047 [arXiv:1110.2865 [astro-ph.CO]].
  * [66] K. Schmitz, “New Sensitivity Curves for Gravitational-Wave Experiments” [Data set], Zenodo, [http://doi.org/10.5281/zenodo.3689582].
  * [67] K. Schmitz, “New Sensitivity Curves for Gravitational-Wave Signals from Cosmological Phase Transitions,” JHEP 01 (2021), 097 [arXiv:2002.04615 [hep-ph]].
  * [68] S. Kuroyanagi, K. Nakayama and J. Yokoyama, “Prospects of determination of reheating temperature after inflation by DECIGO,” PTEP 2015 (2015) no.1, 013E02 [arXiv:1410.6618 [astro-ph.CO]].
  * [69] T. Asaka, S. Blanchet and M. Shaposhnikov, “The nuMSM, dark matter and neutrino masses,” Phys. Lett. B 631 (2005) 151 [hep-ph/0503065].
  * [70] T. Asaka and M. Shaposhnikov, “The nuMSM, dark matter and baryon asymmetry of the Universe,” Phys. Lett. B 620 (2005) 17 [hep-ph/0505013].
  * [71] W. Buchmuller, V. Domcke and K. Schmitz, “Spontaneous B-L Breaking as the Origin of the Hot Early Universe,” Nucl. Phys. B 862 (2012), 587-632 [arXiv:1202.6679 [hep-ph]].
  * [72] W. Buchmuller, V. Domcke, K. Kamada and K. Schmitz, “A Minimal Supersymmetric Model of Particle Physics and the Early Universe,” [arXiv:1309.7788 [hep-ph]].
  * [73] W. Buchmuller, V. Domcke, H. Murayama and K. Schmitz, “Probing the scale of grand unification with gravitational waves,” Phys. Lett. B 809 (2020), 135764 [arXiv:1912.03695 [hep-ph]].
  * [74] Y. Ema, K. Hamaguchi, T. Moroi and K. Nakayama, “Flaxion: a minimal extension to solve puzzles in the standard model,” JHEP 01 (2017), 096 [arXiv:1612.05492 [hep-ph]].
  * [75] P. Campeti, E. Komatsu, D. Poletti and C. Baccigalupi, “Measuring the spectrum of primordial gravitational waves with CMB, PTA and Laser Interferometers,” [arXiv:2007.04241 [astro-ph.CO]].
  * [76] T. Hiramatsu, M. Kawasaki, T. Sekiguchi, M. Yamaguchi and J. Yokoyama, “Improved estimation of radiated axions from cosmological axionic strings,” Phys. Rev. D 83 (2011), 123531 [arXiv:1012.5502 [hep-ph]].
  * [77] T. Hiramatsu, M. Kawasaki, K. Saikawa and T. Sekiguchi, “Production of dark matter axions from collapse of string-wall systems,” Phys. Rev. D 85 (2012), 105020 [arXiv:1202.5851 [hep-ph]].
  * [78] M. Kawasaki, K. Saikawa and T. Sekiguchi, “Axion dark matter from topological defects,” Phys. Rev. D 91 (2015) no.6, 065014 [arXiv:1412.0789 [hep-ph]].
  * [79] V. B. Klaer and G. D. Moore, “How to simulate global cosmic strings with large string tension,” JCAP 10 (2017), 043 [arXiv:1707.05566 [hep-ph]].
  * [80] M. Gorghetto, E. Hardy and G. Villadoro, “Axions from Strings: the Attractive Solution,” JHEP 07 (2018), 151 [arXiv:1806.04677 [hep-ph]].
  * [81] A. Vaquero, J. Redondo and J. Stadler, “Early seeds of axion miniclusters,” JCAP 04 (2019), 012 [arXiv:1809.09241 [astro-ph.CO]].
  * [82] M. Buschmann, J. W. Foster and B. R. Safdi, “Early-Universe Simulations of the Cosmological Axion,” Phys. Rev. Lett. 124 (2020) no.16, 161103 [arXiv:1906.00967 [astro-ph.CO]].
  * [83] M. Hindmarsh, J. Lizarraga, A. Lopez-Eiguren and J. Urrestilla, “Scaling Density of Axion Strings,” Phys. Rev. Lett. 124 (2020) no.2, 021301 [arXiv:1908.03522 [astro-ph.CO]].
  * [84] A. Drew and E. P. S. Shellard, “Radiation from Global Topological Strings using Adaptive Mesh Refinement: Methodology and Massless Modes,” [arXiv:1910.01718 [astro-ph.CO]].
  * [85] N. Du et al. [ADMX], “A Search for Invisible Axion Dark Matter with the Axion Dark Matter Experiment,” Phys. Rev. Lett. 120 (2018) no.15, 151301 [arXiv:1804.05750 [hep-ex]].
  * [86] W. Chung, “CULTASK, The Coldest Axion Experiment at CAPP/IBS in Korea,” PoS CORFU2015 (2016), 047
  * [87] L. Zhong et al. [HAYSTAC], “Results from phase 1 of the HAYSTAC microwave cavity axion experiment,” Phys. Rev. D 97 (2018) no.9, 092001 [arXiv:1803.03690 [hep-ex]].
  * [88] A. Caldwell et al. [MADMAX Working Group], “Dielectric Haloscopes: A New Way to Detect Axion Dark Matter,” Phys. Rev. Lett. 118 (2017) no.9, 091801 [arXiv:1611.05865 [physics.ins-det]].
  * [89] P. Brun et al. [MADMAX], “A new experimental approach to probe QCD axion dark matter in the mass range above 40 $\mu$eV,” Eur. Phys. J. C 79 (2019) no.3, 186 [arXiv:1901.07401 [physics.ins-det]].
  * [90] A. A. Melcón et al., “Axion Searches with Microwave Filters: the RADES project,” JCAP 05 (2018), 040 [arXiv:1803.01243 [hep-ex]].
  * [91] M. Quiros, “Finite temperature field theory and phase transitions,” [arXiv:hep-ph/9901312 [hep-ph]].
  * [92] M. E. Carrington, “The Effective potential at finite temperature in the Standard Model,” Phys. Rev. D 45 (1992), 2933-2944.