Induced Gravitational Waves from Ultra Slow-Roll Inflation
and Pulsar Timing Arrays Observations

Hassan Firouzjahi<EMAIL_ADDRESS>Alireza Talebian<EMAIL_ADDRESS>
1School of Astronomy, Institute for Research in Fundamental Sciences (IPM)
P. O. Box 19395-5531, Tehran, Iran

###### Abstract

The stochastic gravitational wave background (SGWB) detected recently by the
pulsar timing arrays (PTAs) observations may have cosmological origins. In
this work we consider a model of single field inflation containing an
intermediate phase of ultra slow-roll. Fixing the amplitude of the peak of
curvature perturbations by the PBHs bounds we calculate the gravitational
waves (GWs) induced from the curvature perturbations enhanced during USR. The
spectrum of the induced GWs depends on the sharpness of the transition from
the USR phase to the final attractor phase as well as to the duration of the
USR period. While the model can accommodate the current PTAs data but it has
non-trivial predictions for the induced GWs on higher frequency ranges which
can be tested by future observations.

##  1 Introduction

There are indications of detection of stochastic gravitational waves
background (SGWB) from recent pulsar timing arrays (PTAs) around the frequency
range $\sim 10~{}\mathrm{nHz}$ as reported in NANOGrav [1], Parkers PTA [2],
European PTA [3] and the China PTA [4]. These signals may have cosmological
origins as well as astrophysical interpretations. A natural astrophysical
interpretation of the observed SGWB is the superpositions of gravitational
waves (GWs) signals from the merging of binary supermassive black holes [1].
On the other hand, if the observed signal has cosmological origins, this can
open a new window into observing the primordial universe and to probe physics
beyond the Standard Model of particle physics. Among possible cosmological
interpretations of the SGWB are the GWs induced from the enhanced scalar
perturbations on small scales generated during inflation, first order
cosmological phase transitions [5, 6, 7], domain walls or cosmic strings [8,
9, 10, 11, 12], see [13, 14] for further review. It should be noted that the
previous NANOGrav 12.5 year data [15] also indicated traces of SGWB with a
flat spectrum in a narrow range of $\mathrm{nHz}$ frequency which initiated
interests to look for the origins of this signal [16, 17, 18, 19].

Scalar induced gravitational waves (SIGW) from the enhanced scalar
perturbations on small scale during inflation [20, 21, 22, 23, 24, 25, 26, 27,
28, 29] is a mechanism which can explain the observed SGWBs [13, 14]. In this
mechanism, the GWs are sourced at the second order in perturbation theory via
their interaction with the scalar sector generated during inflation.
Typically, this setup requires that the amplitude of scalar perturbations to
grow by about seven orders of magnitude compared to the observed CMB scale.
Consequently, this mechanism can yield to primordial black holes (PBHs)
formation which may comprise all or parts of dark matter energy density [30,
31, 32, 33, 34]. Naively speaking, suppose $k_{\rm p}$ is the peak scale of
perturbations responsible for PBH formation. Then PBH with a given mass have
an induced GW counterpart with a peak frequency at around

$\displaystyle\dfrac{f_{\rm p}}{\rm 10nHz}\simeq 1.54\times
10^{-7}\bigg{(}\dfrac{k_{\rm p}}{{\rm Mpc}^{-1}}\bigg{)}\,.$ (1.1)

For example, in the case of PTAs signal, the amplified mode $k_{\rm p}$ leaves
the horizon during inflation about $\sim 16$ e-fold after the CMB scale with
$k_{\rm CMB}=0.05~{}{\rm Mpc}^{-1}$ has crossed the horizon.

The setup of ultra slow-roll (USR) inflation has been employed as a model in
which the primordial power spectrum can be enhanced to induce large SIGWs and
PBHs [35, 36, 37, 38, 39, 40, 41, 42, 43], for a review see [33, 34]. The USR
setup is a single field model of inflation in which the potential is flat [44,
45, 46] and the inflaton velocity falls off exponentially so the curvature
perturbations grow on superhorizon scales [47]. Since the curvature
perturbations grow on superhorizon scales the USR setup violates the Maldacena
non-Gaussianity consistency condition [48, 49] in single field models [47, 50,
51, 52, 53, 54, 55, 56, 57]. Originally, it was shown in [47] that the
amplitude of local-type non-Gaussianity in USR setup is $f_{NL}=\frac{5}{2}$.
This question was further examined in [58] in which it was shown that the
amplitude of $f_{NL}$ crucially depends on the sharpness of the transition
from the USR phase to the final slow-roll (SR) phase. In an extreme sharp
transition from the USR phase to the SR phase, which was assumed in [47],
$f_{NL}$ reaches its maximum value $\frac{5}{2}$. However, if the transition
to the final stage is mild then the curvature perturbations evolve after the
USR phase before it reaches to its final attractor value. Correspondingly, in
a mild transition, the amplitude of $f_{NL}$ is washed out in the subsequent
evolution and it ends up with a value at the order of the slow-roll
parameters.

Another important point is the issue of loop corrections in this setup. This
question was studied in various recent works [59, 60, 61, 62, 63, 64, 65, 66,
67, 68, 69, 70, 71]. Originally, it was argued in [59], see also [60], that
loop corrections from small scale modes which leave the horizon during the USR
phase can induce large corrections on CMB scale modes. This was criticized in
[61, 62] arguing, among other things, that for a mild transition the loop
corrections will be less significant and the standard PBHs formation within
the single field USR scenario is applicable. This question was studied in some
details in [65] with emphasis on the effects of the sharpness of the
transition from the intermediate USR phase to the final attractor SR phase. It
was shown in [65] that for an arbitrarily sharp transition the one-loop
corrections can be very large, in line with the results advocated in [59, 60].
However, it was speculated in [65] that for a mild transition, the dangerous
one-loop corrections are washed out during the subsequent evolution of the
mode function after the USR phase. This conclusion was further examined in
[66] confirming the physical expectations of [61, 62]. In summary, in order
for the one-loop corrections on CMB scale modes to be harmless one needs a
mild enough transition from the USR phase to the final attractor phase.

In this paper we employ the USR setup as a working mechanism to generate large
SIGW as a possible explanations for the observed SGWB in the NANOGrav data
[1]. For various recent works on SIGWs as an explanation of the the PTAs data
see [73, 74, 75, 76, 72, 77, 78, 79, 80, 81, 82].

##  2 The Setup

The setup which we use to enhance the primordial power spectrum to induce
large GWs at the range of scales observed by the PTA observations contains an
intermediate phase of USR in single field inflation. We have a three-stage
model of inflation in which the large CMB scale leaves the horizon at the
early stage of inflation. The first stage of inflation proceeds say in about
16 e-folds or so. Then the USR phase takes over in which the potential is very
flat and the curvature perturbations experience a growth on superhorizon
scales. In order for the curvature power spectrum to be under perturbative
controls the USR phase has to be terminated followed by a final phase of SR
inflation. A realistic setup requires that the transition from the first SR to
the USR stage and from the USR phase to the final SR phase to be smooth.
However, in order to follow the dynamics analytically, we consider an
idealized situation in which the transition from the SR to USR and then to
final SR phase are instantaneous. Assuming that USR phase is extended during
the time interval $t_{i}\leq t\leq t_{e}$, we assume that the transitions at
the starting point $t=t_{i}$ and at the final point $t=t_{e}$ to be
instantaneous. While the transition to the final SR phase is instantaneous,
but it will take time for the system to relax to its final attractor phase. We
control this relaxation period by a sharpness parameter which plays important
role in our analysis below. It it important that the instantaneous gluing of
the solutions should not be confused with the sharpness of the transition to
the final attractor solution.

With the above discussions in mind, let us elaborate on the dynamics of our
setup. During the first stage of inflation, $t<t_{i}$, the system follows an
attractor phase and the dynamics of the inflaton field $\phi$ is given by the
usual SR dynamics. The Hubble expansion rate $H\equiv\frac{\dot{a}}{a}$ is
nearly constant in which $a(t)$ is the FLRW scale factor. The small evolution
of $H$ during the first stage of inflation is measured by the first SR
parameter $\epsilon\equiv-\frac{\dot{H}}{H^{2}}$ which is nearly constant and
small. During the USR phase the dynamics of the system is given by

$\displaystyle\ddot{\phi}+3H\dot{\phi}=0\,,\quad\quad 3M_{P}^{2}H^{2}\simeq
V_{0}\,,$ (2.1)

where $M_{P}$ is the reduced Planck mass.

As the potential is flat during the USR phase, $H$ approaches a fixed value
and from the field equation we obtain $\dot{\phi}\propto a(t)^{-3}$.
Correspondingly, the slow-roll parameter $\epsilon$ falls off exponentially
during the USR phase as well, $\epsilon\propto a(t)^{-6}$. On the other hand,
the second slow-roll parameter
$\eta\equiv\frac{\dot{\epsilon}}{H\epsilon}\simeq-6$ which is the hallmark of
the USR phase.

It is convenient to work with the number of e-fold $N$ as the clock,
$dN=H(t)dt$. We choose the convention that the USR phase starts at $N=0$ so
for the first stage of inflation, $N<0$. In particular, the CMB scale modes
leave the horizon at around $N\sim-16$. The duration of the USR phase is
denoted by $\Delta N$ which is a free parameter of our setup. Going to
conformal time, $d\tau=dt/a(t)$, the USR phase is extended during
$\tau_{i}\leq\tau\leq\tau_{e}$ and the duration of USR phase is given by
$\Delta
N=\ln\big(\frac{\tau_{i}}{\tau_{e}}\big{missing})=\ln\big(\frac{k_{e}}{k_{i}}\big{missing})$
in which $k_{i}(k_{e})$ represents the mode which leaves the horizon at the
start (end) of USR phase. The slow-roll parameter at the end of USR phase
$\epsilon_{e}$ is related to its value at the start of USR phase
$\epsilon_{i}$ via $\epsilon_{e}=\epsilon_{i}e^{-6\Delta N}$.

As explained above, we assume the USR phase is followed by a SR phase.
Therefore, we need to investigate the evolution of the slow-roll parameters
$\epsilon(\tau)$ and $\eta(\tau)$ after the USR phase. This was studied in
more details in [58], see also [83, 84] for similar studies. Let us denote the
final SR parameters with their attractor values by $\epsilon_{\mathrm{V}}$ and
$\eta_{\mathrm{V}}$ which are expressed in terms of the first and the second
derivatives of the potential in the final SR phase. To simplify the analysis,
here, as in [65], we assume that the potential in the final stage is such that
$\epsilon_{\mathrm{V}}\gg|\eta_{\mathrm{V}}|$ though this assumption can be
relaxed with no important changes in the results. A key parameter in our setup
is the sharpness parameter $h$ which controls how quickly the system reaches
to its final attractor limit. Following [58], we define $h$ as follows

$\displaystyle
h\equiv\frac{6\sqrt{2\epsilon_{\mathrm{V}}}}{\dot{\phi}(t_{e})}=-6\sqrt{\frac{\epsilon_{V}}{\epsilon_{{e}}}}\,.$
(2.2)

With this definition, the slow-roll parameters $\epsilon(\tau)$ and
$\eta(\tau)$ after the USR transition are given by

$\displaystyle\epsilon(\tau)=\epsilon_{e}\Big{[}\frac{h}{6}-(1+\frac{h}{6})\big{(}\frac{\tau}{\tau_{e}}\big{)}^{3}\Big{]}^{2}\,\quad\quad(\tau>\tau_{e})\,,$
(2.3)

and

$\displaystyle\eta(\tau)=-\frac{6(6+h)}{(6+h)-h\big{(}\frac{\tau_{e}}{\tau}\big{)}^{3}}\,\quad\quad(\tau>\tau_{e}).$
(2.4)

As discussed previously, the above results are obtained in the limit of an
instant transition from the USR phase to the final SR phase. Even in this
limit it will take some time for the system to reach to its attractor phase
which is measured by the sharpness parameter $h$. In the limit
$h\rightarrow-\infty$, the system reaches its final attractor phase
immediately after $\tau_{e}$ in which the mode functions become frozen. On the
other hand, for other values of $h$ the system keeps evolving after $\tau_{e}$
until $\epsilon(\tau)$ approaches its final attractor value
$\epsilon_{\mathrm{V}}$. A particular case of transition is when $h=-6$ in
which $\epsilon(\tau)$ is frozen to its value at the end of USR,
$\epsilon(\tau)=\epsilon_{e}$ with $\eta(\tau)=0$ for $\tau>\tau_{e}$. This
limit was mostly studied in recent literature concerning the loop correction
such as in [59]. In the following analysis, as in [65], we consider a general
value of $h$ including the spacial case of $h=-6$. Another important point is
that the larger the value of $|h|$ is the larger $\epsilon_{V}$ is compared to
$\epsilon_{e}$. Correspondingly, the final power spectrum scales somewhat
inversely with $|h|$. As a result, a larger (smaller) value of $|h|$ yields to
a smaller (larger) final power spectrum.

We work with the comoving curvature perturbation ${\cal{R}}$ which in
spatially flat gauge is related to inflaton perturbation via
${\cal{R}}\equiv-\frac{H}{\dot{\phi}}\delta\phi$. Going to Fourier space, we
define the quantum mode function in terms of the annihilation and creation
operators as usual via

$\displaystyle{\cal{R}}(t,{\bf x})=\int\frac{d^{3}k}{(2\pi)^{3}}e^{i{\bf
k}\cdot{\bf x}}\Big{(}{\cal{R}}_{k}(t)a_{\bf k}+{\cal{R}}^{*}_{k}(t)a_{-\bf
k}^{\dagger}\Big{)}\,,$ (2.5)

in which $a_{\bf k}$ and $a_{\bf k}^{\dagger}$ satisfy the usual commutation
relation associated to the annihilation and creation operators, $[a_{\bf
k},a_{{\bf k}^{\prime}}^{\dagger}]=\delta({\bf k}-{\bf k}^{\prime})$.

Starting with the Bunch-Davies (Minkowski) initial condition and imposing the
continuity of ${\cal{R}}$ and $\dot{\cal{R}}$ at the transition points
$\tau=\tau_{{i}}$ and $\tau=\tau_{{e}}$, the mode function at each stage of
$SR\rightarrow USR\rightarrow SR$ can be obtained [65]. The outgoing curvature
perturbation ${\cal{R}}^{(3)}_{k}(t)$ in the final USR phase (third phase) is
given by [65],

$\displaystyle{\cal{R}}^{(3)}_{k}=\frac{H}{M_{P}\sqrt{4\epsilon(\tau)k^{3}}}\Big{[}\alpha^{(3)}_{k}(1+ik\tau)e^{-ik\tau}+\beta^{(3)}_{k}(1-ik\tau)e^{ik\tau}\Big{]}\,,$
(2.6)

with $\epsilon(\tau)$ given by Eq. (2.3) and the coefficients
$(\alpha^{(3)}_{k},\beta^{(3)}_{k})$ are as follow,

$\displaystyle\alpha^{(3)}_{k}=\frac{1}{8k^{6}\tau_{i}^{3}\tau_{e}^{3}}\Big{[}3h(1-ik\tau_{e})^{2}(1+ik\tau_{i})^{2}e^{2ik(\tau_{e}-\tau_{i})}-i(2k^{3}\tau_{i}^{3}+3ik^{2}\tau_{i}^{2}+3i)(4ik^{3}\tau_{e}^{3}-hk^{2}\tau_{e}^{2}-h)\Big{]},$

and

$\displaystyle\beta^{(3)}_{k}=\frac{1}{8k^{6}\tau_{i}^{3}\tau_{e}^{3}}\Big{[}3(1+ik\tau_{i})^{2}(h+hk^{2}\tau_{e}^{2}+4ik^{3}\tau_{e}^{3})e^{-2ik\tau_{i}}+ih(1+ik\tau_{e})^{2}(3i+3ik^{2}\tau_{i}^{2}+2k^{3}\tau_{i}^{3})e^{-2ik\tau_{e}}\Big{]}.$

The physical quantities are calculated at the end of inflation $\tau=0$ when
the system has reached to its attractor phase with
$\epsilon(\tau)\rightarrow\epsilon_{V}$. The curvature perturbation power
spectrum ${\cal{P}}_{\cal{R}}$ from Eq. (2.6) is obtained to be

$\displaystyle{\cal{P}}_{\cal{R}}(k,\tau=0)=\frac{H^{2}}{8M_{P}^{2}\pi^{2}\epsilon_{\mathrm{V}}}\big{|}\alpha^{(3)}_{k}+\beta^{(3)}_{k}\big{|}^{2}\,.$
(2.7)

The behaviour of ${\cal{P}}_{\cal{R}}(k,\tau=0)$ are plotted in Fig. 1. There
are a number of common features as can be seen in this figure. First, we have
the plateau on large scales associated to the modes which leave the horizon
long before the USR phase starts. The amplitude of the power spectrum for
these perturbations is fixed by the COBE normalization on
$\mathrm{k}_{\text{\tiny CMB}}=0.05\,\mathrm{Mpc}^{-1}$ with
${\cal{P}}_{\cal{R}}\simeq 2.1\times 10^{-9}$. Second, prior to the USR phase
there is a dip in power spectrum followed by a universal scaling
${\cal{P}}\propto k^{4}$. Third, there are oscillations superimposed on the
USR plateau after the maximum. Fourth, there is a plateau for the modes which
leave the horizon at the final stage of inflation. As discussed previously,
the larger is the value of $|h|$, the smaller is the final power spectrum
which we will demonstrate bellow.

Figure 1:  The plot of ${\cal{P}}_{\cal{R}}$ vs. $k$ for various values of $h$
with $h=-0.1$ (red), $h=-1$ (orange), $h=-6$ (green) and $h=-12$ (blue). Left:
The values of $(h,\Delta N)$ are fixed such that the peak of
${\cal{P}}_{\cal{R}}$ does not violate the PBHs bounds shown by the black
(top) curve. All four curves very nearly share the same values of the USR peak
and the dip prior to the USR phase. Right: We have fixed $\Delta N=1.3$. As
$\Delta N$ is fixed for all curves, the power is controlled by $h$ such that
the larger $|h|$ is the smaller the final plateau is. In addition, the
position of the dip moves towards the right by increasing $|h|$.

Let us look at the parameter dependence of the the power spectrum given in Eq.
(2.7). Two important parameters are the sharpness of the transition $h$ and
the duration of the USR phase $\Delta N$. In addition, we have the energy
scale of inflation $H$ and the final slow-roll parameter $\epsilon_{V}$. As
can be seen from Eq. (2.7) the latter two parameters appear in a combination
which is fixed by the overall COBE normalization at CMB scales. Another
independent variable may be considered to be the starting time of USR phase,
$\tau_{i}$. However, in order to obtain the enhancement in power for PTAs
observations, we need the starting time of USR to be when the mode of interest
which leaves the horizon has the nano-Hertz frequency. This requires the
starting time of USR phase compared to the CMB scales to be separated by about
16 e-folds. Finally, the spectral index $n_{s}$ is fixed by its best fit value
from Planck observation, i.e. $n_{s}\simeq 0.96$ [85, 86]. In summary, at the
end we have two main independent parameters: the sharpness parameter $h$ and
the duration of USR $\Delta N$ which we will vary in our analysis below.

A crucial point is that models with the intermediate USR phase can generate
significant PBHs which are constrained by observations. Therefore, in order to
meet the various bounds on PBHs formation, we need to impose an additional
bound on $(h,\Delta N)$ parameter space. These PBH constraints leave only one
of them free which we take to be $h$. A view of the power spectrum for various
values of $h$ and the bound from PBHs are shown in Fig. 1. Schematically, we
see that the PBHs bound is roughly translated into
${\cal{P}}_{\cal{R}}<10^{-2}$.

More precisely, in order to consider the PBHs bound on the curvature power
spectrum, we need to know about the probability distribution function (PDF) of
the primordial perturbations. In a common approach [32], the mass fraction
parameter $\beta$ is related to statistics of ${\cal R}$ as [87]

$\displaystyle\beta\simeq\int_{{\cal R}_{c}}^{\infty}~{}f_{\cal
R}(x)~{}\differential x\simeq\dfrac{1}{2}{\rm Erfc}\left(\dfrac{{\cal
R}_{c}}{\sqrt{2{\cal P}_{\cal R}}}\right)$ (2.8)

where $f_{\cal R}$ is PDF of $\cal R$ and ${\cal R}_{c}\sim{\cal O}(1)$ [88,
33]. The second estimation comes from the fact that we can consider a Gaussian
PDF for $\cal R$ with zero mean and variance at the order of power spectrum.

After PBH production, it is crucial to determine the fraction of PBH abundance
in dark matter density at the present epoch. It is roughly given by [32]

$f_{\text{PBH}}(M_{\rm PBH})\simeq 2.7\times 10^{8}\Big{(}\dfrac{M_{\rm
PBH}}{M_{\odot}}\Big{)}^{-\frac{1}{2}}\beta(M_{\rm PBH})\,,$ (2.9)

where $M_{\odot}$ and $M_{\rm PBH}$ are the solar mass and the PBH mass
respectively. Assuming an instant reheating at the end of inflation [33], PBH
mass can be estimated by

$\displaystyle\dfrac{M_{\rm PBH}}{M_{\odot}}\simeq 30\,\bigg{(}\dfrac{k_{\rm
p}}{3.2\times 10^{5}~{}{\rm Mpc}^{-1}}\bigg{)}^{-2}\,,$ (2.10)

where $k_{\rm p}$ is the peak scale of perturbations responsible for PBH
formation, i.e. the wave-number associated to the maximum of the power
spectrum. In Fig. 2, we have illustrated the mass function $f_{\rm PBH}$ for
various values of $k_{\rm p}$. PBHs with a given mass have an induced GW
counterpart with a peak frequency at around Eq. (1.1). Therefore, we have a
one-to-one correspondence between the frequency and the PBH mass.

Figure 2:  Fraction $f_{\rm PBH}$ as a function of the mass of the formed PBHs
in unit of solar mass for USR models in which the power spectrum has a peak
around $k_{\rm p}$. The observational constraints are taken from Refs. [95,
96, 97]. The parameters $(h,\Delta N)$ are fixed such that the curve for each
$k_{\mathrm{p}}$ reaches the maximum value of power allowed by the PBHs bound.
For example, for $k_{\mathrm{p}}=2.3\times 10^{7}~{}{\rm Mpc}^{-1}$ (green
curve), we have $(h,\Delta N)=(-0.1,\,1.45),(-1,\,2.17),(-6,\,2.59)$ and
$(-12,\,2.68)$. These values of $(h,\Delta N)$ are used to generate the
results in Figs. 3 and 4.

Now let us look at the various limits of the power spectrum Eq. (2.7). We have
two dimensionless numbers, $x\equiv-k\tau_{i}$ and $e^{-\Delta N}$. First
consider the limit $e^{-\Delta N}\ll x$ so we expand Eq. (2.7) to leading
order in $e^{-\Delta N}$, obtaining

$\displaystyle{\cal{P}}_{\cal{R}}(k,\tau=0)\simeq$
$\displaystyle\frac{e^{6\Delta N}}{2}{\cal P}_{\text{\tiny
CMB}}\big{(}\frac{h-6}{h}\big{)}^{2}$ (2.11)
$\displaystyle\times\Big{[}2x^{6}+9x^{4}+18x^{2}+9+(21x^{4}-9)\cos(2x)+(6x^{5}-24x^{3}-18x)\sin(2x)\Big{]}\,,$

in which ${\cal P}_{\text{\tiny CMB}}$ is the CMB scale power spectrum given
by

$\displaystyle{\cal P}_{\text{\tiny
CMB}}=\frac{H^{2}}{8\pi^{2}M_{P}^{2}\epsilon_{i}}\,.$ (2.12)

From Eq. (2.11) we see that ${\cal{P}}_{\cal{R}}\propto e^{6\Delta N}$ which
is the hallmark of USR inflation for the modes which leave the horizon during
the USR phase. Second, we see that
${\cal{P}}_{\cal{R}}\propto\big{(}\frac{h-6}{h}\big{)}^{2}$. This is clearly
seen in Fig. 1 as cases with higher value of $|h|$ have lower power in the
final plateau. The physical reason is as follows. Models with larger $|h|$
reach the final attractor phase more quickly. For this to happen,
$\epsilon(\tau)$ should assume its final value
$\epsilon_{\mathrm{V}}>\epsilon_{e}$ quickly as well. This means that the mode
becomes frozen quickly after the USR phase but with a final amplitude
${\cal{P}}_{\cal{R}}(k,\tau=0)<{\cal{P}}_{\cal{R}}(k,\tau_{e})$.

To understand the scaling behaviour of the power spectrum prior to USR phase
and close to the USR peak, let us consider the $x\ll 1$ limit of the
expression (2.11), obtaining

$\displaystyle{\cal{P}}_{\cal{R}}(k,\tau=0)\simeq\frac{2}{25}e^{6\Delta
N}{\cal P}_{\text{\tiny
CMB}}\big{(}\frac{h-6}{h}\big{)}^{2}\,(k\tau_{i})^{4}\,.$ (2.13)

It shows that the power spectrum scales like ${\cal{P}}_{\cal{R}}(k)\propto
k^{4}$ prior to and after the USR phase starts, a phenomenon which was
observed previously in [89, 90, 91, 92, 93] as well.

As we see in Fig. 1, there is a dip in power spectrum prior to USR phase where
the above mentioned $k^{4}$ scaling starts. To understand the nature of this
dip, note that the expression (2.11) is obtained assuming that $e^{-\Delta
N}\ll x$. However, this limit is violated for very long modes which become
superhorizon much earlier than the USR phase starts. In particular, the CMB
scale modes belong to this limit. Considering the $x\ll e^{-\Delta N}$ limit
of the power spectrum we obtain

$\displaystyle{\cal{P}}_{\cal{R}}(k,\tau=0)\simeq{\cal P}_{\text{\tiny
CMB}}\Big{(}1-\frac{4}{5}\frac{h-6}{h}\,(k\tau_{i})^{2}\Big{)}\,,\quad\quad(k\tau_{i}\rightarrow
0)\,.$ (2.14)

The position of the dip $k=k_{\mathrm{d}}$ is where the two expressions (2.14)
and (2.13) become comparable, yielding to the approximate value (see also
[93])

$\displaystyle
k_{\mathrm{d}}\tau_{i}\simeq\sqrt{\frac{5h}{4(h-6)}}\,e^{-\frac{3}{2}\Delta
N}\,.$ (2.15)

From the above formula we see that for a fixed value of $\Delta N$, as $|h|$
increase the value of $k_{\mathrm{d}}$ increases as well, i.e. the dip moves
towards the right, as seen in the right panel of Fig. 1.

As we mentioned previously the USR model can generate non-Gaussianities.
However, the magnitude of $f_{NL}$ depends on $k$ as well. For the mode which
leaves the horizon during the early stage of inflation and prior to USR phase,
then the Maldacena consistency condition does hold and for these modes
$f_{NL}$ is basically very small. On the other hand, for the modes which leave
the horizon during the USR phase, i.e. for $k_{i}<k<k_{e}$, the consistency
condition is violated. The final value of $f_{NL}$ for these modes crucially
depends on the parameter $h$. This was studied in details in [58] and [65] in
which it is shown that up to slow-roll corrections,

$\displaystyle f_{NL}=\frac{5h^{2}}{2(h-6)^{2}}\,.$ (2.16)

For an infinitely sharp transition with $h\rightarrow-\infty$ in which the
system assumes its final attractor phase immediately after the USR transition,
we obtain the maximum value $f_{NL}=\frac{5}{2}$. However, lowering $|h|$
reduces $f_{NL}$ accordingly. For example, for the standard scenario in which
$h=-6$ as studied in [59] one obtains $f_{NL}=\frac{5}{8}\simeq 0.63$. For
milder transitions with $|h|\lesssim 1$, from Eq. (2.16) we typically obtain
$f_{NL}\ll 1$. For example for $h=-1$ and $h=-0.1$ which we will study below,
we obtain $f_{NL}\simeq 0.051$ and $f_{NL}\simeq 0.0007$ respectively.
Therefore, to very good approximation one can employ the Gaussian bound on
PBH’s formation. To be more specific, to neglect the non-Gaussianity effects
in PBH formation, we need that $f_{NL}{\cal{P}}_{\cal{R}}\ll 1$ [94]. In our
model with the maximum value $f_{NL}=\frac{5}{2}$, we can easily satisfy the
Gaussianity approximation if ${\cal{P}}_{\cal{R}}$ is small. In our analysis,
as can be seen in Fig. 1, the PBHs bound typically require that
${\cal{P}}_{\cal{R}}\lesssim 10^{-2}$ so we easily meet the condition
$f_{NL}{\cal{P}}_{\cal{R}}\ll 1$ for all practical values of $h$.

As mentioned in Introduction section, the loop correction is an open question
in this setup [59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71]. The
effects of the sharpness of the transition on the loop corrections were
studied in [65]. It was shown that for an arbitrarily sharp transition the
one-loop corrections become very large. More specifically, it was shown in
[65] that for $|h|\gg 1$, the one-loop corrections scale linearly with $h$,
yielding to a large loop correction in line with the results of [59, 60].
However, it was shown in [66] that for a mild transition the one-loop
corrections to CMB scale modes are slow-roll suppressed and are harmless. In
our current study, in order to trust the Gaussian approximation for PBHs
formation and to neglect the large loop corrections, one requires $|h|\lesssim
1$ which we shall assume in our analysis. However, in order to compare the
predictions of the setup with both sharp and mild transitions, we also present
the results for SIGWs for the cases of sharp transition as well. In our
numerical plots below, the examples of sharp transitions correspond to the
cases $h=-6$ and $h=-12$. As another perspective to the question of loop
corrections, it was argued in [71] that by considering the surface terms the
leading loop corrections will disappear. This was studied for the cubic
Hamiltonian in [71] while the contribution associated to quartic Hamiltonian
and the surface terms is an open question.

##  3 SIGWs and Constraints from PTAs Observations

The curvature perturbations, which have been generated during inflation, can
re-enter the horizon during the radiation-dominated (RD) era in which the
metric (in conformal Newtonian gauge) reads

$\displaystyle\differential s^{2}=$
$\displaystyle-a^{2}\left[(1+2\Phi)\differential\tau^{2}+\Big{(}(1-2\Psi)\delta_{ij}+\dfrac{1}{2}h_{ij}\Big{)}\differential
x^{i}\differential x^{j}\right].$ (3.1)

Here $\tau$ is the conformal time during RD era, $\Phi$ and $\Psi$ are the
Bardeen potentials and $h_{ij}$ is the transverse-traceless tensor
perturbations. Using the Einstein’s field equations, the evolution of Fourier
modes of $h_{ij}$, denoted by $h^{\lambda}_{\bf{k}}$, are given by

$\displaystyle{h^{\lambda}_{\bf{k}}}^{\prime\prime}\\!(\eta)+2\mathcal{H}{h^{\lambda}_{\bf{k}}}^{\prime}\\!(\eta)+k^{2}h^{\lambda}_{\bf{k}}(\eta)=4S^{\lambda}_{\bf{k}}(\eta),$
(3.2)

where $\lambda$ represents the two polarizations. The primes denote the
derivative with respect to the conformal time, $\mathcal{H}=a^{\prime}/a$ is
the conformal Hubble parameter, and the source term $S^{\lambda}_{\bf k}$ is
transverse and traceless which is second order in scalar perturbations, given
by

$\displaystyle S^{\lambda}_{\bf{k}}=$
$\displaystyle\,\int\\!\frac{\differential^{3}q}{(2\pi)^{3}}\,\varepsilon^{\lambda}_{ij}(\hat{\bf
k})~{}q^{i}q^{j}\bigg{[}2\Phi_{\bf{q}}\Phi_{\bf{k-q}}+\left(\mathcal{H}^{-1}\Phi^{\prime}_{\bf{q}}+\Phi_{\bf{q}}\right)\left(\mathcal{H}^{-1}\Phi^{\prime}_{\bf{k-q}}+\Phi_{\bf{k-q}}\right)\bigg{]}\,,$
(3.3)

where $\varepsilon^{\lambda}_{ij}$ is the polarization tensor. Note that here
we have neglected the vector perturbations and the anisotropic stress
($\Phi\simeq\Psi$).

In Fourier space, the Bardeen potential is related to ${\cal R}_{\bf k}$
through transfer function ${\cal T}({\bf k}\tau)$ as

$\displaystyle\Phi_{\bf k}=\dfrac{2}{3}{\cal T}({\bf k}\tau){\cal R}_{\bf
k}\,.$ (3.4)

The transfer function encodes the linear evolution of the Newtonian potential
after horizon reentry which has a oscillatory behaviour. Solving the equation
of motion (3.2), taking the late-time limit during a RD era
($\tau\rightarrow\infty$ at the matter-radiation equality), the power spectrum
of tensor fluctuations is given by [98]

$\displaystyle{\cal P}_{h}(\tau,k)=4\int_{0}^{\infty}\differential
v\int_{|1-v|}^{|1+v|}\differential
u~{}~{}\mathcal{K}\left(u,v,k\tau\right)~{}\mathcal{P}_{\mathcal{\cal
R}}\left(ku\right)\mathcal{P}_{\mathcal{\cal R}}\left(kv\right)\,,$ (3.5)

For further details about the integration kernel $\mathcal{K}$ and how to
perform the integrations see [98].

The produced GW contributes to the total energy of Universe and dilutes like
radiation. Taking into account the following matter-dominated and dark-energy-
dominated eras, the current value of $\Omega_{\rm GW}$, the fraction of the GW
energy density per logarithmic wavelength, is obtained to be

$\displaystyle\Omega_{\rm GW}h_{0}^{2}=\Omega_{\rm
r}h_{0}^{2}~{}\left(\dfrac{g_{*}}{g_{*,e}}\right)^{1/3}\Omega_{{\rm
GW},e}(f)\,.$ (3.6)

Here $\Omega_{\rm r}$ is the present-day abundance of radiation, $g_{*}$ is
the number of relativistic degrees of freedom in energy density, and the
subscripts $e$ denotes the time of emission. Note that $\Omega_{\rm
r}h_{0}^{2}\simeq 4.2\times 10^{-5}$ with $h_{0}=H_{0}/100\,{\rm km}\,{\rm
s}^{-1}\,{\rm Mpc}^{-1}$. Here $f=c\,k/(2\pi)$ is the frequency of the GW
which has appeared due to the $k$-dependence of the spectrum of curvature
perturbations (2.7) during inflation.

We have used the curvature perturbations power spectrum (2.7) generating
during the USR phase and calculated the convolution integral (3.5) numerically
to find the current fractional energy density of GWs (3.6). The results are
shown in Fig. 3 for various values of $h$ and $\Delta N$ in nano-Hertz bound.
The results have been presented alongside the posteriors of an Hellings-Downs
(HD) correlated free spectral reconstruction of the NANOGrav signal [13]. The
values of $(h,\Delta N)$ are fixed such that the peak of ${\cal P}_{\cal R}$
does not violate the PBHs bounds. As seen, the spectra of our model follow the
HD pattern expected for a stochastic gravitational wave background.
Interestingly, we see that within the observable window the predictions of all
models, mild ($|h|\lesssim 1$) or sharp ($|h|\gg 1$), follow the same pattern
and are not significantly different from each other. However, outside the
NANOGrav observed window ($f>10^{2}\,\mathrm{nHz}$) the curves deviate from
each other noticeably. This pattern is similar to the plots of power spectrum
of curvature perturbations presented in Fig. 1. The reason is that all curves
are obtained after imposing the PBHs bounds. However, the starting time of the
USR and the value of the peak of the USR plateau are very similar for all
curves as seen in Fig. 1. This is the reason why all curves, sharp or mild,
follow close trajectories on the observable window.

However, crucial to our setup is that outside the NANOGrav window, the curves
have distinct predictions for SIGWs on frequencies much higher than $\sim
10^{2}\,\mathrm{nHz}$. More specifically, the final tail of the power spectrum
scales somewhat inversely with the sharpness parameter $h$ such that milder
(sharper) transitions have higher (lower) tails. In Fig. 4 we we have shown
the SIGWs spectra for a larger frequency range. In this figure, the quantity
$\Omega_{\rm GW}h^{2}$ was plotted against the frequency together with the
sensitivity of the various existing and forthcoming GW experiments such as
LISA, SKA, BBO, DECIGO etc. As seen, the tails of SIGW spectrums for different
$(h,\Delta N)$ can fall into the sensitivity bounds of these observations. It
means that different values of $(h,\Delta N)$ are distinguishable from each
other in future GWs observations.

Figure 3:  The prediction of the model compared to the NANOGrav data [13] for
a power spectrum with the peak scale $k_{\mathrm{p}}\simeq 10^{7}~{}{\rm
Mpc}^{-1}$ but for various values of $h$. The value of $\Delta N$ is fixed by
the PBHs fraction allowed in Fig. 2. More specifically, $(h,\Delta
N)=(-0.1,\,1.45),(-1,\,2.17),(-6,\,2.59)$ and $(-12,\,2.68)$.  Figure 4:  The
same as in Fig. 3 but for an extended range of the frequency. As we see,
depending on the value of $h$, the tail of SIGWs differ noticeably for each
model which can be tested in future observations. The higher is the value of
$|h|$ the lower is the value of the SIGW tail.

##  4 Summary and Discussions

The stochastic gravitational wave background detected by various PTAs
observations can open a new window to the dynamics of early universe. In
particular, this signal can be generated by GWs induced by scalar
perturbations at second order in perturbation theory. The SIGWs can be used as
a tool to distinguish various inflationary scenarios. A key criteria is that
the models which are employed to explain the SGWB observed by the PTAs
observations should not generate too much of PBHs which are constrained in
various frequency ranges.

In this work we have considered a single field model of inflation containing
an intermediate USR phase. This setup has been used extensively in the past to
generate PBHs and for the induced GWs studies. We have paid particular
attention to the sharpness parameter of the model which play significant roles
in loop corrections and for the amplitude of non-Gaussianity. In order to be
away from the dangerous one-loop corrections we require a mild transition with
$|h|\lesssim 1$. This is also the limit where the amplitude of non-Gaussianity
is small and one can employ the Gaussian predictions of the PBHs formation.

We have calculated the spectrum of SIGWs and compared it to the NANOGrave
results. The predictions of the model are consistent with the observed data.
However, a careful data analysis is required to contrast the predictions of
the model with the PTA datas and to put constraints on the model parameters.

While our setup can qualitatively explain the origin of the NANOGrave
observations, but it has specific predictions for the spectrum in higher
frequency ranges. Our model predicts that a setup with a mild (sharp)
transition has a higher (lower) tail of IGWs once they are fit to the current
NANOGrave data. As a result, the predictions of our model for the amplitude of
induced GWs can be tested in future GWs observation which may put constraints
on the model parameter or to rule it out.

Acknowledgments: We are grateful to Anotonio Riotto, Sina Hooshangi, Seyed Ali
Hosseini Mansoori, and Ali Akbar Abolhasani for useful discussions and
comments on the draft. We are partially supported by the “Saramadan”
Federation of Iran. A. T. would like to thank University of Rwanda, EAIFR, and
ICTP for their kind hospitalities during the 17th international workshop on
the “Dark Side of the Universe” where this work was in its final stages.

## References

  * [1] G. Agazie et al. [NANOGrav], Astrophys. J. Lett. 951, no.1, L8 (2023) doi:10.3847/2041-8213/acdac6 [arXiv:2306.16213 [astro-ph.HE]].
  * [2] D. J. Reardon et al, Astrophys. J. Lett. 951, no.1, L6 (2023) doi:10.3847/2041-8213/acdd02 [arXiv:2306.16215 [astro-ph.HE]].
  * [3] J. Antoniadis et al, [arXiv:2306.16214 [astro-ph.HE]].
  * [4] H. Xu et al, doi:10.1088/1674-4527/acdfa5 [arXiv:2306.16216 [astro-ph.HE]].
  * [5] A. Kosowsky, M. S. Turner and R. Watkins, Phys. Rev. Lett. 69, 2026-2029 (1992).
  * [6] M. Kamionkowski, A. Kosowsky and M. S. Turner, Phys. Rev. D 49, 2837-2851 (1994).
  * [7] C. Caprini, R. Durrer and G. Servant, Phys. Rev. D 77, 124015 (2008).
  * [8] M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, Phys. Rev. Lett. 112, 041301 (2014).
  * [9] T. W. B. Kibble, J. Phys. A 9, 1387-1398 (1976) doi:10.1088/0305-4470/9/8/029
  * [10] A. Vilenkin, Phys. Lett. B 107, 47-50 (1981) doi:10.1016/0370-2693(81)91144-8
  * [11] R. R. Caldwell and B. Allen, Phys. Rev. D 45, 3447-3468 (1992) doi:10.1103/PhysRevD.45.3447
  * [12] A. Vilenkin, Phys. Rev. D 23, 852-857 (1981) doi:10.1103/PhysRevD.23.852
  * [13] A. Afzal et al. [NANOGrav], Astrophys. J. Lett. 951, no.1, L11 (2023).
  * [14] J. Antoniadis et al, [arXiv:2306.16227 [astro-ph.CO]].
  * [15] Z. Arzoumanian et al. [NANOGrav], Astrophys. J. Lett. 905, no.2, L34 (2020).
  * [16] S. Vagnozzi, Mon. Not. Roy. Astron. Soc. 502, no.1, L11-L15 (2021).
  * [17] Z. Zhou, J. Jiang, Y. F. Cai, M. Sasaki and S. Pi, Phys. Rev. D 102, no.10, 103527 (2020)
  * [18] G. Domènech and S. Pi, Sci. China Phys. Mech. Astron. 65, no.3, 230411 (2022)
  * [19] V. Dandoy, V. Domcke and F. Rompineve, [arXiv:2302.07901 [astro-ph.CO]].
  * [20] K. N. Ananda, C. Clarkson and D. Wands, Phys. Rev. D 75, 123518 (2007).
  * [21] D. Baumann, P. J. Steinhardt, K. Takahashi and K. Ichiki, Phys. Rev. D 76, 084019 (2007).
  * [22] E. Bugaev and P. Klimai, Phys. Rev. D 81, 023517 (2010).
  * [23] H. Assadullahi and D. Wands, Phys. Rev. D 79, 083511 (2009).
  * [24] L. Alabidi, K. Kohri, M. Sasaki and Y. Sendouda, JCAP 09, 017 (2012).
  * [25] R. g. Cai, S. Pi and M. Sasaki, Phys. Rev. Lett. 122, no.20, 201101 (2019).
  * [26] S. Pi and M. Sasaki, JCAP 09, 037 (2020).
  * [27] S. Balaji, G. Domenech and J. Silk, JCAP 09, 016 (2022).
  * [28] A. Talebian, S. A. Hosseini Mansoori and H. Firouzjahi, Astrophys. J. 948, no.1, 48 (2023)
  * [29] G. Domènech, Universe 7, no.11, 398 (2021).
  * [30] B. Carr, F. Kuhnel and M. Sandstad, Phys. Rev. D 94, no.8, 083504 (2016).
  * [31] B. Carr and F. Kuhnel, Ann. Rev. Nucl. Part. Sci. 70, 355-394 (2020).
  * [32] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, Class. Quant. Grav. 35, no.6, 063001 (2018).
  * [33] O. Özsoy and G. Tasinato, [arXiv:2301.03600 [astro-ph.CO]].
  * [34] C. T. Byrnes and P. S. Cole, [arXiv:2112.05716 [astro-ph.CO]].
  * [35] P. Ivanov, P. Naselsky and I. Novikov, Phys. Rev. D 50, 7173-7178 (1994).
  * [36] J. Garcia-Bellido and E. Ruiz Morales, Phys. Dark Univ. 18, 47-54 (2017).
  * [37] M. Biagetti, G. Franciolini, A. Kehagias and A. Riotto, JCAP 07, 032 (2018).
  * [38] H. V. Ragavendra, P. Saha, L. Sriramkumar and J. Silk, Phys. Rev. D 103, no.8, 083510 (2021).
  * [39] H. Di and Y. Gong, JCAP 07, 007 (2018).
  * [40] J. Liu, Z. K. Guo and R. G. Cai, Phys. Rev. D 101, no.8, 083535 (2020)
  * [41] S. Hooshangi, A. Talebian, M. H. Namjoo and H. Firouzjahi, Phys. Rev. D 105, no.8, 083525 (2022)
  * [42] S. Hooshangi, M. H. Namjoo and M. Noorbala, [arXiv:2305.19257 [astro-ph.CO]].
  * [43] A. Ghoshal, A. Moursy and Q. Shafi, [arXiv:2306.04002 [hep-ph]].
  * [44] W. H. Kinney, Phys. Rev. D 72, 023515 (2005) [gr-qc/0503017].
  * [45] M. J. P. Morse and W. H. Kinney, Phys. Rev. D 97, no.12, 123519 (2018).
  * [46] W. C. Lin, M. J. P. Morse and W. H. Kinney, JCAP 09, 063 (2019).
  * [47] M. H. Namjoo, H. Firouzjahi and M. Sasaki, Europhys. Lett. 101, 39001 (2013).
  * [48] J. M. Maldacena, JHEP 0305, 013 (2003) [astro-ph/0210603].
  * [49] P. Creminelli and M. Zaldarriaga, JCAP 10, 006 (2004).
  * [50] J. Martin, H. Motohashi and T. Suyama, Phys. Rev. D 87, no.2, 023514 (2013).
  * [51] X. Chen, H. Firouzjahi, M. H. Namjoo and M. Sasaki, Europhys. Lett. 102, 59001 (2013).
  * [52] X. Chen, H. Firouzjahi, E. Komatsu, M. H. Namjoo and M. Sasaki, JCAP 1312, 039 (2013).
  * [53] M. Akhshik, H. Firouzjahi and S. Jazayeri, JCAP 12, 027 (2015).
  * [54] S. Mooij and G. A. Palma, JCAP 11, 025 (2015).
  * [55] R. Bravo, S. Mooij, G. A. Palma and B. Pradenas, JCAP 05, 024 (2018).
  * [56] B. Finelli, G. Goon, E. Pajer and L. Santoni, Phys. Rev. D 97, no.6, 063531 (2018).
  * [57] S. Pi and M. Sasaki, Phys. Rev. Lett. 131, no.1, 011002 (2023).
  * [58] Y. F. Cai, X. Chen, M. H. Namjoo, M. Sasaki, D. G. Wang and Z. Wang, JCAP 05, 012 (2018).
  * [59] J. Kristiano and J. Yokoyama, arXiv2211.03395hep-th.
  * [60] J. Kristiano and J. Yokoyama, [arXiv:2303.00341 [hep-th]].
  * [61] A. Riotto, [arXiv:2301.00599 [astro-ph.CO]].
  * [62] A. Riotto, [arXiv:2303.01727 [astro-ph.CO]].
  * [63] S. Choudhury, S. Panda and M. Sami, [arXiv:2302.05655 [astro-ph.CO]].
  * [64] S. Choudhury, S. Panda and M. Sami, [arXiv:2303.06066 [astro-ph.CO]].
  * [65] H. Firouzjahi, [arXiv:2303.12025 [astro-ph.CO]].
  * [66] H. Firouzjahi and A. Riotto, arXiv2304.07801astro-ph.CO.
  * [67] H. Motohashi and Y. Tada, [arXiv:2303.16035 [astro-ph.CO]].
  * [68] G. Franciolini, A. Iovino, Junior., M. Taoso and A. Urbano, arXiv2305.03491astro-ph.CO.
  * [69] G. Tasinato, arXiv2305.11568hep-th.
  * [70] H. Firouzjahi, [arXiv:2305.01527 [astro-ph.CO]].
  * [71] J. Fumagalli, [arXiv:2305.19263 [astro-ph.CO]].
  * [72] S. Wang, Z. C. Zhao, J. P. Li and Q. H. Zhu, [arXiv:2307.00572 [astro-ph.CO]].
  * [73] S. Vagnozzi, [arXiv:2306.16912 [astro-ph.CO]].
  * [74] G. Franciolini, A. Iovino, Junior., V. Vaskonen and H. Veermae, [arXiv:2306.17149 [astro-ph.CO]].
  * [75] Y. F. Cai, X. C. He, X. Ma, S. F. Yan and G. W. Yuan, [arXiv:2306.17822 [gr-qc]].
  * [76] K. Inomata, K. Kohri and T. Terada, [arXiv:2306.17834 [astro-ph.CO]].
  * [77] L. Liu, Z. C. Chen and Q. G. Huang, [arXiv:2307.01102 [astro-ph.CO]].
  * [78] Z. Yi, Q. Gao, Y. Gong, Y. Wang and F. Zhang, [arXiv:2307.02467 [gr-qc]].
  * [79] D. G. Figueroa, M. Pieroni, A. Ricciardone and P. Simakachorn, [arXiv:2307.02399 [astro-ph.CO]].
  * [80] B. M. Gu, F. W. Shu and K. Yang, [arXiv:2307.00510 [astro-ph.CO]].
  * [81] R. Ebadi, S. Kumar, A. McCune, H. Tai and L. T. Wang, [arXiv:2307.01248 [astro-ph.CO]].
  * [82] E. Madge, E. Morgante, C. P. Ibáñez, N. Ramberg and S. Schenk, [arXiv:2306.14856 [hep-ph]].
  * [83] Y. F. Cai, X. H. Ma, M. Sasaki, D. G. Wang and Z. Zhou, Phys. Lett. B 834, 137461 (2022).
  * [84] Y. F. Cai, X. H. Ma, M. Sasaki, D. G. Wang and Z. Zhou, JCAP 12, 034 (2022).
  * [85] N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)].
  * [86] Y. Akrami et al. [Planck], Astron. Astrophys. 641, A10 (2020).
  * [87] J. Garcia-Bellido, M. Peloso and C. Unal, JCAP 12, 031 (2016).
  * [88] I. Musco, V. De Luca, G. Franciolini and A. Riotto, Phys. Rev. D 103, no.6, 063538 (2021)
  * [89] C. T. Byrnes, P. S. Cole and S. P. Patil, JCAP 06, 028 (2019), [arXiv:1811.11158 [astro-ph.CO]].
  * [90] P. S. Cole, A. D. Gow, C. T. Byrnes and S. P. Patil, [arXiv:2204.07573 [astro-ph.CO]].
  * [91] P. Carrilho, K. A. Malik and D. J. Mulryne, Phys. Rev. D 100, no.10, 103529 (2019), [arXiv:1907.05237 [astro-ph.CO]].
  * [92] O. Özsoy and G. Tasinato, Phys. Rev. D 105, no.2, 023524 (2022).
  * [93] S. Pi and J. Wang, JCAP 06, 018 (2023).
  * [94] S. Young and C. T. Byrnes, JCAP 08, 052 (2013).
  * [95] A. M. Green and B. J. Kavanagh, J. Phys. G 48, no.4, 043001 (2021)
  * [96] Kavanagh, B. J. 2019, bradkav/PBHbounds: Release version, 1.0, Zenodo, doi: 10.5281/zenodo.3538999.
  * [97] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, Rept. Prog. Phys. 84, no.11, 116902 (2021).
  * [98] K. Kohri and T. Terada, Phys. Rev. D 97, no.12, 123532 (2018).


# Inverse magnetic catalysis effect and current quark mass effect on mass
spectra and Mott transitions of pions under external magnetic field

Luyang Li1 Shijun Mao2<EMAIL_ADDRESS>1 Xi’an University of Posts
and Telecommunications, Xi’an, Shaanxi 710121, China
2 Institute of Theoretical Physics, School of Physics, Xi’an Jiaotong
University, Xi’an, Shaanxi 710049, China

###### Abstract

Mass spectra and Mott transition of pions $(\pi^{0},\ \pi^{\pm})$ at finite
temperature and magnetic field are investigated in a two-flavor NJL model, and
we focus on the inverse magnetic catalysis (IMC) effect and current quark mass
(CQM) effect. Due to the dimension reduction of the constituent quarks, the
pion masses jump at their Mott transitions, which is independent of the IMC
effect and CQM effect. We consider the IMC effect by using a magnetic
dependent coupling constant, which is a monotonic decreasing function of
magnetic field. With IMC effect, the Mott transition temperature of $\pi^{0}$
meson $T_{m}^{0}$ is a monotonic decreasing function of magnetic field. For
charged pions $\pi^{\pm}$, the Mott transition temperature $T_{m}^{+}$ fast
increases in weak magnetic field region and then decreases with magnetic
field, which are accompanied with some oscillations. Comparing with the case
without IMC effect, $T_{m}^{0}$ and $T_{m}^{+}$ are lower when including IMC
effect. CQM effect are considered by varying parameter $m_{0}$ in non-chiral
limit. For $\pi^{0}$ meson, $T_{m}^{0}$ is not a monotonic function of
magnetic field with low $m_{0}$, but it is a monotonic decreasing function
with larger $m_{0}$. In the weak magnetic field region, $T_{m}^{0}$ is higher
for larger $m_{0}$, but in the strong magnetic field region, it is lower for
larger $m_{0}$. For $\pi^{+}$ meson, $T^{+}_{m}$ is only quantitatively
modifies by current quark mass effect, and it becomes higher with larger
$m_{0}$.

##  I Introduction

The study of hadron properties in QCD medium is important for our
understanding of strong interaction matter, due to its close relation to QCD
phase structure and relativistic heavy ion collision. For instance, the chiral
symmetry breaking leads to the rich meson spectra, and the mass shift of
hadrons will enhance or reduce their thermal production in relativistic heavy
ion collisions k1 ; k2 ; k3 .

It is widely believed that the strongest magnetic field in nature may be
generated in the initial stage of relativistic heavy ion collisions. The
initial magnitude of the field can reach $eB\sim(1-100)m_{\pi}^{2}$ in
collisions at the Relativistic Heavy Ion Collider and the Large Hadron
Collider b0 ; b1 ; b2 ; b3 ; b4 , where $e$ is the electron charge and
$m_{\pi}$ the pion mass in vacuum. In recent years, the study of magnetic
field effects on the hadrons attract much attention. As the Goldstone bosons
of the chiral (isospin) symmetry breaking, the properties of neutral (charged)
pions at finite magnetic field, temperature and density are widely
investigated c1 ; c3 ; hadron1 ; hadron2 ; qm1 ; qm2 ; sigma1 ; sigma2 ;
sigma3 ; sigma4 ; l1 ; l2 ; l3 ; l4 ; lqcd5 ; lqcd6 ; ding2008.00493 ;
ding2022 ; njl2 ; meson ; mfir ; ritus5 ; ritus6 ; mao1 ; mao11 ; mao2 ; wang
; coppola ; phi ; liuhao3 ; he ; Restrepo ; maocharge ; maopion ;
yulang2010.05716 ; q1 ; q2 ; q3 ; q4 ; huangamm1 ; huangamm2 ; q5 ; q6 ; q7 ;
q8 ; q9 ; q10 ; su3meson4 ; Scoccola1 ; Scoccola2 ; meijie3 .

The LQCD simulations performed with physical pion mass observe the inverse
magnetic catalysis (IMC) phenomenon lattice1 ; lattice11 ; lattice2 ; lattice3
; lattice4 ; lattice5 ; lattice6 ; lattice7 . Namely, the pseudo-critical
temperature $T_{pc}$ of chiral symmetry restoration and the quark mass near
$T_{pc}$ drop down with increasing magnetic field. On analytical side, many
scenarios are proposed to understand this IMC phenomenon, but the physical
mechanism is not clear fukushima ; mao ; hadron2 ; bf1 ; bf12 ; bf13 ; bf2 ;
bf3 ; bf4 ; bf5 ; bf51 ; bf52 ; bf8 ; bf9 ; bf11 ; db1 ; db2 ; db3 ; db5 ; db6
; pnjl1 ; pnjl2 ; pnjl3 ; pnjl4 ; pqm ; ferr1 ; ferr2 ; huangamm2 ; meijie ;
imcreview ; ammmao . Besides, how does the IMC effect influence the pion
properties under external magnetic field? Previous studies focus on the case
with vanishing temperature and density, and report that IMC effect leads to a
lower mass for pions under external magnetic field meson ; su3meson4 ; mao11 .

The current quark mass determines the explicit breaking of chiral symmetry. In
LQCD simulations, it is considered through different pion mass l4 ; lattice5 ;
lattice3 . At vanishing temperature and finite magnetic field, the normalized
neutral pion mass $m_{\pi^{0}}(eB)/m_{\pi^{0}}(eB=0)$ increases with current
quark mass, but the normalized charged pion mass
$m_{\pi^{\pm}}(eB)/m_{\pi^{\pm}}(eB=0)$ decreases with current quark mass l4 .
In effective models, the current quark mass (CQM) effect on pion properties
has not yet been studied.

In this paper, we will investigate the IMC effect and CQM effect on mass
spectra and Mott transition of pions at finite temperature and magnetic field.
Here we make use of a Pauli-Villars regularized Nambu–Jona-Lasinio (NJL)
model, which describes remarkably well the static properties of light mesons
njl1 ; njl2 ; njl3 ; njl4 ; njl5 ; zhuang . In our calculations, the IMC
effect is introduced by a magnetic field dependent coupling, and the CQM
effect is considered by tuning the quark mass parameter.

The rest of paper is organized as follows. Sec.II introduces our theoretical
framework for meson spectra and the Mott transition in a Pauli-Villars
regularized NJL model. The numerical results and discussions are presented in
Sec.III, which focus on inverse magnetic catalysis effect in Sec.IIIA and
current quark mass effect in Sec.IIIB. Finally, we give the summary in Sec.IV.

##  II Framework

The two-flavor NJL model is defined through the Lagrangian density in terms of
quark fields $\psi$ njl1 ; njl2 ; njl3 ; njl4 ; njl5 ; zhuang

${\cal
L}=\bar{\psi}\left(i\gamma_{\nu}D^{\nu}-m_{0}\right)\psi+G\left[\left(\bar{\psi}\psi\right)^{2}+\left(\bar{\psi}i\gamma_{5}{\vec{\tau}}\psi\right)^{2}\right].$
(1)

Here the covariant derivative $D_{\nu}=\partial_{\nu}+iQA_{\nu}$ couples
quarks with electric charge $Q=diag(Q_{u},Q_{d})=diag(2e/3,-e/3)$ to the
external magnetic field ${\bf B}=(0,0,B)$ in $z$-direction through the
potential $A_{\nu}=(0,0,Bx_{1},0)$. $m_{0}$ is the current quark mass, which
determines the explicit breaking of chiral symmetry. $G$ is the coupling
constant in scalar and pseudo-scalar channels, which determines the
spontaneously breaking of chiral symmetry and isospin symmetry.

In NJL model, mesons are constructed through quark bubble summations in the
frame of random phase approximation njl2 ; njl3 ; njl4 ; njl5 ; zhuang ,

${\cal D}_{M}(x,z)=2G\delta(x-z)+\int d^{4}y\ 2G\Pi_{M}(x,y){\cal
D}_{M}(y,z),$ (2)

where ${\cal D}_{M}(x,y)$ represents the meson propagator from $x$ to $y$ in
coordinate space, and the corresponding meson polarization function is the
quark bubble,

$\Pi_{M}(x,y)=i{\text{Tr}}\left[\Gamma_{M}^{\ast}S(x,y)\Gamma_{M}S(y,x)\right]$
(3)

with the meson vertex

$\Gamma_{M}=\left\\{\begin{array}[]{ll}1&M=\sigma\\\
i\tau_{+}\gamma_{5}&M=\pi_{+}\\\ i\tau_{-}\gamma_{5}&M=\pi_{-}\\\
i\tau_{3}\gamma_{5}&M=\pi_{0}\
,\end{array}\right.\Gamma_{M}^{*}=\left\\{\begin{array}[]{ll}1&M=\sigma\\\
i\tau_{-}\gamma_{5}&M=\pi_{+}\\\ i\tau_{+}\gamma_{5}&M=\pi_{-}\\\
i\tau_{3}\gamma_{5}&M=\pi_{0}\ ,\end{array}\right.$

the quark propagator matrix in flavor space $S=diag(S_{u},\ S_{d})$, and the
trace in spin, color and flavor spaces.

According to the Goldstone’s theorem, the pseudo-Goldstone mode of chiral
(isospin) symmetry breaking under external magnetic field is the neutral pion
$\pi^{0}$ (charged pions $\pi^{\pm}$) gold1 ; gold2 . The charged pions are no
longer the pseudo-Goldstone modes since their direct interaction with magnetic
field.

####  II.0.1 neutral pion $\pi^{0}$

The neutral pion $\pi^{0}$ is affected by external magnetic field only through
the pair of charged constituent quarks, and its propagator in momentum space
can be derived as

${\cal D}_{\pi^{0}}(k)=\frac{2G}{1-2G\Pi_{\pi^{0}}(k)},$ (4)

with the polarization function $\Pi_{\pi^{0}}(k)$ and conserved momentum
$k=(k_{0},{\bf k})$ of $\pi^{0}$ meson under external magnetic field.

The meson pole mass $m_{\pi^{0}}$ is defined as the pole of the propagator at
zero momentum ${\bf k}={\bf 0}$,

$1-2G\Pi_{\pi^{0}}(\omega^{2}=m^{2}_{\pi^{0}},{\bf k}^{2}=0)=0.$ (5)

At nonzero magnetic field, the three-dimensional quark momentum integration in
the polarization function $\Pi_{\pi^{0}}$ becomes a one-dimensional momentum
integration and a summation over the discrete Landau levels. The polarization
function can be simplified

$\displaystyle\Pi_{\pi^{0}}(\omega^{2},0)=J_{1}(m_{q})+\omega^{2}J_{2}(\omega^{2})$
(6)

and

$\displaystyle J_{1}(m_{q})$ $\displaystyle=$ $\displaystyle
3\sum_{f,n}\alpha_{n}\frac{|Q_{f}B|}{2\pi}\int\frac{dp_{3}}{2\pi}\frac{1-2F(E_{f})}{E_{f}},$
$\displaystyle J_{2}(\omega^{2})$ $\displaystyle=$ $\displaystyle
3\sum_{f,n}\alpha_{n}\frac{|Q_{f}B|}{2\pi}\int\frac{dp_{3}}{2\pi}{{1-2F(E_{f})}\over
E_{f}(4E_{f}^{2}-w^{2})},$

with the summation over all flavors and Landau energy levels, spin factor
$\alpha_{n}=2-\delta_{n0}$, quark energy
$E_{f}=\sqrt{p^{2}_{3}+2n|Q_{f}B|+m_{q}^{2}}$, and Fermi-Dirac distribution
function $F(x)=\left(e^{x/T}+1\right)^{-1}$.

The (dynamical) quark mass $m_{q}$ is determined by the gap equation,

$\displaystyle 1-2GJ_{1}(m_{q})$ $\displaystyle=$
$\displaystyle\frac{m_{0}}{m_{q}}.$ (7)

During the chiral symmetry restoration, the quark mass decreases, and the
$\pi^{0}$ mass increases, as guaranteed by the Goldstone’s theorem gold1 ;
gold2 . When the $\pi^{0}$ mass is beyond the threshold

$m_{\pi^{0}}=2m_{q},$ (8)

the decay channel $\pi^{0}\rightarrow q{\bar{q}}$ opens, which defines the
${\pi^{0}}$ Mott transition mott1 ; zhuang ; mott2 ; mott3 .

From the explicit expression of $\Pi_{\pi^{0}}$ in Eq(6), the factor
$1/(4E_{f}^{2}-\omega^{2})$ in the integrated function of $J_{2}(\omega^{2})$
becomes $(1/4)/(p_{3}^{2}+2n|Q_{f}B|)$ at $\omega=2m_{q}$. When we do the
integration over $p_{3}$, the $p_{3}^{2}$ in the denominator leads to the
infrared divergence at the lowest Landau level $n=0$. Therefore,
$m_{\pi^{0}}=2m_{q}$ is not a solution of the pole equation, and there must be
a mass jump for $\pi^{0}$ meson at the Mott transition. This mass jump is a
direct result of the quark dimension reduction mao2 ; maocharge ; q2 ; q3 ;
yulang2010.05716 , and independent of the parameters and regularization
schemes of NJL model. When the magnetic field disappears, there is no more
quark dimension reduction, the integration $\int d^{3}{\bf
p}/(4E_{f}^{2}-\omega^{2})\sim\int dp$ becomes finite at $\omega=2m_{q}$, and
there is no more such a mass jump. It should be mentioned that in chiral
limit, no such mass jump happens for ${\pi^{0}}$ meson even under external
magnetic field, which has been analytically proved in our previous work mao1 .

####  II.0.2 charged pions $\pi^{\pm}$

When constructing charged mesons through quark bubble summations, we should
take into account of the interaction between charged mesons and magnetic
fields. The charged pions $\pi^{\pm}$ with zero spin are Hermite conjugation
to each other, they have the same mass at finite temperature and magnetic
field.

The $\pi^{+}$ meson propagator $D_{\pi^{+}}$ can be expressed in terms of the
polarization function $\Pi_{\pi^{+}}$ maocharge ; mao11 ; maopion ,

$\displaystyle
D_{\pi^{+}}({\bar{k}})=\frac{2G}{1-2G\Pi_{\pi^{+}}({\bar{k}})},$ (9)

where ${\bar{k}}=(k_{0},0,-\sqrt{(2l+1)eB},k_{3})$ is the conserved Ritus
momentum of $\pi^{+}$ meson under magnetic fields.

The meson pole mass $m_{\pi^{+}}$ is defined through the pole of the
propagator at zero momentum $(l=0,\ k_{3}=0)$,

$\displaystyle 1-2G\Pi_{\pi^{+}}(k_{0}=m_{\pi^{+}})=0,$ (10)

and

$\displaystyle\Pi_{\pi^{+}}(k_{0})$ $\displaystyle=$ $\displaystyle
J_{1}(m_{q})+J_{3}(k_{0}),$ (11) $\displaystyle J_{3}(k_{0})$ $\displaystyle=$
$\displaystyle\sum_{n,n^{\prime}}\int\frac{dp_{3}}{2\pi}\frac{j_{n,n^{\prime}}(k_{0})}{4E_{n}E_{n^{\prime}}}\times$
$\displaystyle\left[\frac{F(-E_{n^{\prime}})-F(E_{n})}{k_{0}+E_{n^{\prime}}+E_{n}}+\frac{F(E_{n^{\prime}})-F(-E_{n})}{k_{0}-E_{n^{\prime}}-E_{n}}\right],$
$\displaystyle j_{n,n^{\prime}}(k_{0})$ $\displaystyle=$
$\displaystyle\left[{(k_{0})^{2}/2}-n^{\prime}|Q_{u}B|-n|Q_{d}B|\right]j^{+}_{n,n^{\prime}}$
(13) $\displaystyle-2\sqrt{n^{\prime}|Q_{u}B|n|Q_{d}B|}\
j^{-}_{n,n^{\prime}},$

with $u$-quark energy
$E_{n^{\prime}}=\sqrt{p^{2}_{3}+2n^{\prime}|Q_{u}B|+m_{q}^{2}}$, $d$-quark
energy $E_{n}=\sqrt{p^{2}_{3}+2n|Q_{d}B|+m_{q}^{2}}$, and summations over
Landau levels of $u$ and $d$ quarks in $J_{3}(k_{0})$.

The quark dimension reduction also leads to infrared ($p_{3}\to 0$)
singularity of the quark bubble $\Pi_{\pi^{+}}(k_{0})$ at some Landau level
and some temperature mao1 ; mao2 , and thus there is no solution of the
corresponding pole equation for the $\pi^{+}$ meson mass $m_{\pi^{+}}$ in this
case. Because the spins of $u$ and ${\bar{d}}$ quarks at the lowest-Landau-
level are always aligned parallel to the external magnetic field, and
$\pi^{+}$ meson has spin zero. The two constituents at the lowest Landau level
($n=n^{\prime}=0$) do not satisfy the pole equation (10). Namely, they can not
form a charged $\pi^{+}$ meson. Considering $|Q_{d}|<|Q_{u}|$, the threshold
for the singularity of $\Pi_{\pi^{+}}$ is located at Landau levels
$n^{\prime}=0$ and $n=1$,

$\displaystyle m_{\pi^{+}}=m_{q}+\sqrt{2{|Q_{d}B|}+m_{q}^{2}}.$ (14)

This defines the Mott transition of $\pi^{+}$ meson, and a mass jump will
happen. There exist other mass jumps located at $n^{\prime}\geq 1,n\geq 0$,
see examples in Fig.6 and Fig.12. All these mass jumps are caused by the quark
dimension reduction under external magnetic field, and do not depend on the
parameters and regularization schemes of NJL model.

##  III results and discussions

Because of the four-fermion interaction, the NJL model is not a renormalizable
theory and needs regularization. In this work, we make use of the gauge
invariant Pauli-Villars regularization scheme njl1 ; njl2 ; njl3 ; njl4 ; njl5
; mao ; maocharge ; mao2 , where the quark momentum runs formally from zero to
infinity. The three parameters in the Pauli-Villars regularized NJL model,
namely the current quark mass $m_{0}=5$ MeV, the coupling constant $G=3.44$
GeV-2 and the Pauli-Villars mass parameter $\Lambda=1127$ MeV are fixed by
fitting the chiral condensate $\langle\bar{\psi}\psi\rangle=-(250\
\text{MeV})^{3}$, pion mass $m_{\pi}=134$ MeV and pion decay constant
$f_{\pi}=93$ MeV in vacuum with $T=\mu=0$ and $eB=0$. In our current
calculations, we consider the situation with finite temperature and magnetic
field and vanishing quark chemical potential $\mu=0$.

###  III.1 Inverse magnetic catalysis effect

From LQCD simulations, the inverse magnetic catalysis phenomenon can be
characterized either by the chiral condensates or the pseudo-critical
temperature of chiral symmetry restoration lattice1 ; lattice11 ; lattice2 ;
lattice3 ; lattice4 ; lattice5 ; lattice6 ; lattice7 . Therefore, to include
the inverse magnetic catalysis effect in the NJL model, one approach is to fit
the LQCD results of chiral condensates geb1 ; meson ; geb3 ; geb4 , and
another approach is to fit the LQCD result of pseudo-critical temperature bf8
; bf9 ; geb1 ; su3meson4 ; mao11 .

In our calculations, following Refs bf8 ; bf9 ; geb1 ; su3meson4 ; mao11 , we
use a two-flavor NJL model with a magnetic field dependent coupling $G(eB)$,
which is determined by fitting the LQCD reported decreasing pseudo-critical
temperature of chiral symmetry restoration $T_{pc}(eB)/T_{pc}(eB=0)$ lattice1
. As plotted in Fig.1, the magnetic field dependent coupling $G(eB)/G(eB=0)$
is a monotonic decreasing function of magnetic field, and it reduces $16\%$ at
$eB/m^{2}_{\pi}=30$. In this paper, we fix $m_{\pi}=134$ MeV as the scale of
magnetic field. As we have checked, with our fitted coupling constant $G(eB)$,
the magnetic catalysis phenomena of chiral condensates at low temperature and
the inverse magnetic catalysis phenomena at high temperature can be
reproduced.

Figure 1: Magnetic field dependent coupling $G(eB)$ fitted from LQCD reported
decreasing pseudo-critical temperature of chiral symmetry restoration
$T_{pc}(eB)/T_{pc}(eB=0)$.

####  III.1.1 neutral pion $\pi^{0}$

Figure 2: $\pi^{0}$ meson mass $m_{\pi^{0}}$ as a function of magnetic field
at vanishing temperature, with (red solid line) and without (blue dashed line)
inverse magnetic catalysis effect.

With our fitted $G(eB)$ in Fig.1, we solve the gap equation (7) and pole
equation (5) to obtain the $\pi^{0}$ meson mass at finite temperature and
magnetic field with IMC effect. The results with and without IMC effect are
plotted in red and blue, respectively.

Figure 2 depicts the $\pi^{0}$ meson mass $m_{\pi^{0}}$ as a function of
magnetic field at vanishing temperature with (red solid line) and without
(blue dashed line) inverse magnetic catalysis effect. Because in both cases
the magnetic field enhances the breaking of chiral symmetry in vacuum. As the
pseudo-Goldstone boson, the $\pi^{0}$ meson masses are decreasing functions of
magnetic field with and without IMC effect. We observe a lower value for
$m_{\pi^{0}}$, when including the IMC effect. Similar conclusion is obtained
in Refs meson ; su3meson4 , where inverse magnetic catalysis effect is
introduced into the two-flavor and three-flavor NJL models.

Figure 3: $\pi^{0}$ meson mass $m_{\pi^{0}}$ as a function of temperature
under fixed magnetic field $eB=20m_{\pi}^{2}$, with (red lines) and without
(blue lines) inverse magnetic catalysis effect.

In Fig.3, we plot the $\pi^{0}$ meson mass $m_{\pi^{0}}$ at finite temperature
and fixed magnetic field $eB=20m_{\pi}^{2}$. With (red lines) and without
(blue lines) inverse magnetic catalysis effect, the $\pi^{0}$ mass spectra
have the similar structure. $m_{\pi^{0}}$ increases at low temperature region
and jumps at the Mott transition $T=T^{0}_{m}$. After that, it firstly
decreases and then increases with temperature. The inverse magnetic catalysis
effect leads to some quantitative modifications. For instance, Mott transition
temperature $T^{0}_{m}$ is shifted to a lower value, and $m_{\pi^{0}}$ becomes
lower (higher) at $T<T^{0}_{m}$ ($T>T^{0}_{m}$).

Figure 4: Mott transition temperature of $\pi^{0}$ meson $T^{0}_{m}$ as a
function of magnetic field, with (red solid line) and without (blue dashed
line) inverse magnetic catalysis effect.

Figure 4 shows the Mott transition temperature of $\pi^{0}$ meson $T^{0}_{m}$
as a function of magnetic field, with (red solid line) and without (blue
dashed line) inverse magnetic catalysis effect. Without IMC effect, the
$T^{0}_{m}$ decreases with the magnetic field in weak magnetic field region
$(eB\leq 5m^{2}_{\pi})$, becomes flat in medium magnetic field region
$(5m^{2}_{\pi}\leq eB\leq 16m^{2}_{\pi})$ and decreases in strong magnetic
field region $(eB\geq 16m^{2}_{\pi})$. With IMC effect, the flat structure of
Mott transition temperature disappears and $T^{0}_{m}$ monotonically decreases
with magnetic field. Furthermore, at fixed magnetic field, $T^{0}_{m}$ has a
lower value when including IMC effect.

####  III.1.2 charged pion $\pi^{+}$

With our fitted $G(eB)$ in Fig.1, we solve the gap equation (7) and pole
equation (10) to obtain $\pi^{+}$ meson mass at finite temperature and
magnetic field with IMC effect. The results with and without IMC effect are
plotted in red and blue, respectively.

As shown in Fig.5, with and without IMC effect, $m_{\pi^{+}}$ are increasing
functions of magnetic field, and no decreasing behavior is observed at
vanishing temperature. With IMC effect, $m_{\pi^{+}}$ has a lower value than
without IMC effect, and the deviation becomes larger at stronger magnetic
field. Similar results are obtained in three-flavor NJL model including IMC
effect su3meson4 .

Figure 5: $\pi^{+}$ meson mass $m_{\pi^{+}}$ as a function of magnetic field
at vanishing temperature, with (red solid line) and without (blue dashed line)
inverse magnetic catalysis effect.

In Fig.6, we make comparison of $m_{\pi^{+}}$ as a function of temperature at
fixed magnetic field $eB=20m_{\pi}^{2}$ with (red lines) and without (blue
lines) inverse magnetic catalysis effect. They show similar structure.
$m_{\pi^{+}}$ decreases in the low temperature region. At the Mott transition
$T^{+}_{m}$, $m_{\pi^{+}}$ shows the first jump, and two other mass jumps
happen at $T^{+}_{1}$ and $T^{+}_{2}$. With $T^{+}_{m}<T<T^{+}_{1}$,
$m_{\pi^{+}}$ first decreases and then increases with temperature. With
$T^{+}_{1}<T<T^{+}_{2}$ and $T>T^{+}_{2}$, $m_{\pi^{+}}$ decreases with
temperature. At high enough temperature, $m_{\pi^{+}}$ becomes independent of
temperature and IMC effect. In the low and high temperature regions,
$m_{\pi^{+}}$ with IMC effect is smaller than without IMC effect. The value of
$T^{+}_{m}=124,160$ MeV, $T^{+}_{1}=136,172$ MeV and $T^{+}_{2}=157,187$ MeV
are different for the cases with and without IMC effect, which indicates that
IMC effect lowers down the temperatures of ${\pi^{+}}$ meson mass jumps.

Figure 6: $\pi^{+}$ meson mass $m_{\pi^{+}}$ as a function of temperature at
fixed magnetic field $eB=20m_{\pi}^{2}$, with (red lines) and without (blue
lines) inverse magnetic catalysis effect.

Figure 7 is the Mott transition temperature $T^{+}_{m}$ of $\pi^{+}$ meson as
a function of magnetic field, with (red solid line) and without (blue dashed
line) inverse magnetic catalysis effect. A fast increase of $T^{+}_{m}$ occurs
when turning on external magnetic field, where a peak structure around
$eB\simeq 1m^{2}_{\pi}$ shows up. Without IMC effect, the $T^{+}_{m}$
decreases with magnetic field and then increases, which are associated with
some oscillations. With IMC effect, the $T^{+}_{m}$ decreases as magnetic
field goes up, which is also accompanied with some oscillations. At fixed
magnetic field, $T^{+}_{m}$ has a lower value, when including IMC effect, .

Figure 7: Mott transition temperature of $\pi^{+}$ meson $T^{+}_{m}$ as a
function of magnetic field, with (red solid line) and without (blue dashed
line) inverse magnetic catalysis effect.

###  III.2 Current quark mass effect

In this part, we consider the effect of current quark mass $m_{0}$ on mass
spectra and Mott transition of charged pions. The results are plotted with
black, blue, red, cyan and magenta lines, respectively, corresponding to
$m_{0}=2,\ 5,\ 10,\ 20,\ 50$ MeV. In numerical calculations, we only change
the parameter $m_{0}$ and keep other parameters intact in our NJL model.

####  III.2.1 neutral pion $\pi^{0}$

Figure 8: $\pi^{0}$ meson mass $m_{\pi^{0}}$ (upper panel) and normalized
$\pi^{0}$ meson mass $m_{\pi^{0}}(eB)/m_{\pi^{0}}(eB=0)$ (lower panel) as a
function of magnetic field at vanishing temperature with different current
quark mass $m_{0}$. The results are plotted with black, blue, red, cyan and
magenta lines, respectively, corresponding to $m_{0}=2,\ 5,\ 10,\ 20,\ 50$
MeV.

Figure 8 plots the $\pi^{0}$ meson mass $m_{\pi^{0}}$ and normalized $\pi^{0}$
meson mass $m_{\pi^{0}}(eB)/m_{\pi^{0}}(eB=0)$ as a function of magnetic field
at vanishing temperature with different current quark mass $m_{0}=2,\ 5,\ 10,\
20,\ 50$ MeV. Because the current quark mass determines the explicit breaking
of chiral symmetry. As the pseudo-Goldstone boson, the mass of $\pi^{0}$ meson
will increase with current quark mass when fixing magnetic field and vanishing
temperature. On the other side, magnetic field plays the role of catalysis for
spontaneous breaking of chiral symmetry when fixing current quark mass and
vanishing temperature, and this will lead to a decreasing $m_{\pi^{0}}$. As
shown in Fig.8, with fixed magnetic field, $m_{\pi^{0}}$ becomes larger with
larger $m_{0}$. Moreover, the deviation between $m_{\pi^{0}}$ with different
$m_{0}$ looks independent of the magnetic field. With fixed current quark
mass, $m_{\pi^{0}}$ is a decreasing function of magnetic field. Similar as
$m_{\pi^{0}}$, the normalized $\pi^{0}$ meson mass
$m_{\pi^{0}}(eB)/m_{\pi^{0}}(eB=0)$ increases with the current quark mass when
fixing magnetic field, and decreases with magnetic field when fixing $m_{0}$,
which are consistent with LQCD results l4 .

Figure 9: $\pi^{0}$ meson mass $m_{\pi^{0}}$ (upper panel) and normalized
$\pi^{0}$ meson mass $m_{\pi^{0}}(T)/m_{\pi^{0}}(T=0)$ (lower panel) as a
function of temperature with fixed magnetic field $eB=20m_{\pi}^{2}$ and
different current quark mass $m_{0}$. The results are plotted with black,
blue, red, cyan and magenta lines, respectively, corresponding to $m_{0}=2,\
5,\ 10,\ 20,\ 50$ MeV.

In Fig.9, we depict the $\pi^{0}$ meson mass $m_{\pi^{0}}$ and normalized
$\pi^{0}$ meson mass $m_{\pi^{0}}(T)/m_{\pi^{0}}(T=0)$ as a function of
temperature with fixed magnetic field $eB=20m_{\pi}^{2}$ and different current
quark mass $m_{0}=2,\ 5,\ 10,\ 20,\ 50$ MeV. They show similar structure.
$m_{\pi^{0}}$ slightly increases with temperature, and a mass jump happens at
the Mott transition $T=T^{0}_{m}$. After Mott transition, $m_{\pi^{0}}$ first
decreases and then increases with temperature. However, there exist some
quantitative difference. With larger current quark mass, the $m_{\pi^{0}}$ is
larger in the whole temperature region, and the Mott transition temperature is
lower. At high enough temperature, $m_{\pi^{0}}$ will become degenerate due to
the strong thermal motion of constituent quarks. Different from meson mass
$m_{\pi^{0}}$, the normalized $\pi^{0}$ meson mass
$m_{\pi^{0}}(T)/m_{\pi^{0}}(T=0)$ is larger with smaller current quark mass,
and the deviation between different $m_{0}$ is larger with higher temperature.

Figure 10: Mott transition temperature of $\pi^{0}$ meson $T^{0}_{m}$ as a
function of magnetic field with different current quark mass $m_{0}$. The
results are plotted with black, blue, red, cyan and magenta lines,
respectively, corresponding to $m_{0}=2,\ 5,\ 10,\ 20,\ 50$ MeV.

The Mott transition temperature of $\pi^{0}$ meson $T^{0}_{m}$ is plotted as a
function of magnetic field in Fig.10 with different current quark mass
$m_{0}=2,\ 5,\ 10,\ 20,\ 50$ MeV. For a lower value of current quark mass
$m_{0}=2$ MeV, the Mott transition temperature decreases with weak magnetic
field $(eB\leq 2m^{2}_{\pi})$, slightly increases in medium magnetic field
region $(2m^{2}_{\pi}\leq eB\leq 15m^{2}_{\pi})$ and decreases again in strong
magnetic field region $(eB\geq 15m^{2}_{\pi})$. For $m_{0}=5$ MeV, we obtain a
flat curve for the Mott transition temperature in medium magnetic field region
$(5m^{2}_{\pi}\leq eB\leq 16m^{2}_{\pi})$, and in weak and strong magnetic
field regions, the Mott transition temperature decreases. With the larger
value of current quark mass $m_{0}=10,\ 20,\ 50$ MeV, the Mott transition
temperatures are monotonic decreasing functions of magnetic field. In weak
magnetic field region, $T^{0}_{m}$ is higher for larger $m_{0}$, but in strong
magnetic field region, $T^{0}_{m}$ is lower for larger $m_{0}$.

In the end of this section, we make some comment on chiral limit with
$m_{0}=0$. As discussed in our previous paper mao1 , in chiral limit, the
Goldstone boson $\pi^{0}$ is massless in chiral breaking phase and its mass
continuously increases with temperature in chiral restored phase. The Mott
transition temperature is the same as the critical temperature of chiral
restoration phase transition, and it increases with magnetic field. Moreover,
no mass jump occurs for $\pi^{0}$ meson at the Mott transition with or without
external magnetic field.

####  III.2.2 charged pion $\pi^{+}$

Figure 11: $\pi^{+}$ meson mass $m_{\pi^{+}}$ (upper panel) and normalized
$\pi^{+}$ meson mass $m_{\pi^{+}}(eB)/m_{\pi^{+}}(eB=0)$ (lower panel) as a
function of magnetic field with vanishing temperature and different current
quark mass $m_{0}$. The results are plotted with black, blue, red, cyan and
magenta lines, respectively, corresponding to $m_{0}=2,\ 5,\ 10,\ 20,\ 50$
MeV.

Figure 11 shows $\pi^{+}$ meson mass $m_{\pi^{+}}$ and normalized mass
$m_{\pi^{+}}(eB)/m_{\pi^{+}}(eB=0)$ as a function of magnetic field with
vanishing temperature and different current quark mass $m_{0}=2,\ 5,\ 10,\
20,\ 50$ MeV. $m_{\pi^{+}}$ is an increasing function of magnetic field when
fixing $m_{0}$. With fixed magnetic field, larger value of $m_{0}$ leads to
larger $m_{\pi^{+}}$. In weak magnetic field region, the $m_{0}$ effect is
more obvious than in strong magnetic field region, which is indicated by the
larger difference of $m_{\pi^{+}}$ between different $m_{0}$ cases. However,
with fixed magnetic field, the normalized $\pi^{+}$ meson mass
$m_{\pi^{+}}(eB)/m_{\pi^{+}}(eB=0)$ decreases as $m_{0}$ goes up, which is
consistent with LQCD results l4 .

Figure 12: $\pi^{+}$ meson mass $m_{\pi^{+}}$ as a function of temperature
with fixed magnetic field $eB=20m^{2}_{\pi}$ and different current quark mass
$m_{0}$. The results are plotted with black, blue, red, cyan and magenta
lines, respectively, corresponding to $m_{0}=2,\ 5,\ 10,\ 20,\ 50$ MeV.

$\pi^{+}$ meson mass $m_{\pi^{+}}$ is plotted as a function of temperature
with fixed magnetic field $eB=20m^{2}_{\pi}$ and different current quark mass
$m_{0}=2,\ 5,\ 10,\ 20,\ 50$ MeV in Fig.12. $m_{\pi^{+}}$ has several mass
jumps for all considered current quark mass $m_{0}$, which happen at
temperature $T_{m}^{+}$, $T_{1}^{+}$, $T_{2}^{+}$, $T_{3}^{+}$ successively.
With low temperature $T<T_{m}^{+}$, $m_{\pi^{+}}$ becomes larger with larger
$m_{0}$. With $T_{m}^{+}<T<T_{1}^{+}$ and $T_{1}^{+}<T<T_{2}^{+}$,
$m_{\pi^{+}}$ are also larger with larger $m_{0}$, but the value of
$m_{\pi^{+}}$ is very close for different $m_{0}$. With $T>T_{2}^{+}$,
$m_{\pi^{+}}$ changes nonmonotonically with $m_{0}$. It is clearly shown that
the temperatures with mass jump depend on current quark mass. Within the
considered temperature region $0<T<0.25$ GeV, we observe three jumps in case
of current quark mass $m_{0}=2,\ 5,\ 20$ MeV, four jumps in case of $m_{0}=10$
MeV. For $m_{0}=50$ MeV, the first mass jump occurs at $T=0.262$ GeV, which is
beyond the scope of Fig.12.

Figure 13: Mott transition temperature of $\pi^{+}$ meson $T^{+}_{m}$ as a
function of magnetic field with different current quark mass $m_{0}$. The
results are plotted with black, blue, red, cyan and magenta lines,
respectively, corresponding to $m_{0}=2,\ 5,\ 10,\ 20,\ 50$ MeV.

Figure 13 plots Mott transition temperature of $\pi^{+}$ meson $T^{+}_{m}$ as
a function of magnetic field with different current quark mass $m_{0}=2,\ 5,\
10,\ 20,\ 50$ MeV. $T^{+}_{m}$ shows similar behavior when varying current
quark mass. Accompanied with some oscillations, $T^{+}_{m}$ increases fast in
weak magnetic field region, decreases in medium magnetic field region and
slightly increases in strong magnetic field region. With larger current quark
mass, Mott transition temperature $T^{+}_{m}$ becomes higher. Note that in
case of $m_{0}=50$ MeV, the increasing behavior of $T^{+}_{m}$ appears in very
strong magnetic field, which is beyond the scope of Fig.13.

##  IV summary

Mass spectra and Mott transition of pions $(\pi^{0},\ \pi^{\pm})$ at finite
temperature and magnetic field are investigated in a two-flavor NJL model, and
we focus on the inverse magnetic catalysis effect and current quark mass
effect.

We consider the inverse magnetic catalysis effect by introducing a magnetic
dependent coupling constant into NJL model, which is a monotonic decreasing
function of magnetic field. The mass spectra of pions $(\pi^{0},\ \pi^{\pm})$
at finite temperature and/or magnetic field are not changed qualitatively by
IMC effect. At the Mott transition, the mass jumps of pions happen. Without
IMC effect, the Mott transition temperature of $\pi^{0}$ meson $T^{0}_{m}$
decreases with magnetic field, but shows a flat structure in medium magnetic
field region. With IMC effect, the flat structure of $T^{0}_{m}$ disappears
and $T^{0}_{m}$ is a monotonic decreasing function of magnetic field. For
charged pions $\pi^{\pm}$, the Mott transition temperature $T^{+}_{m}$ is not
a monotonic function of magnetic field. Without IMC effect, it fast increases
in weak magnetic field region, decreases in the medium magnetic field region
and slightly increases in strong magnetic field region, which are accompanied
with some oscillations. When including IMC effect, the increasing behavior of
$T^{+}_{m}$ in strong magnetic field is changed into a decreasing behavior.

The current quark mass effect is studied in non-chiral limit. The masses of
pions $(\pi^{0},\ \pi^{\pm})$ at vanishing temperature increases as the
current quark mass $m_{0}$ goes up. However, the normalized masses
$m_{\pi}(eB)/m_{\pi}(eB=0)$ change differently. For $\pi^{0}$ meson,
$m_{\pi^{0}}(eB)/m_{\pi^{0}}(eB=0)$ increases with $m_{0}$. For $\pi^{+}$
meson, $m_{\pi^{+}}(eB)/m_{\pi^{+}}(eB=0)$ decreases with $m_{0}$. These
properties are consistent with LQCD simulations. At the Mott transition, the
mass jumps of pions happen. The Mott transition temperature of $\pi^{0}$ meson
$T^{0}_{m}$ is qualitatively modifies by current quark mass effect. With a low
value of current quark mass $m_{0}=2$ MeV, $T^{0}_{m}$ decreases with weak
magnetic field, slightly increases in the medium magnetic field region, and
decreases again in the strong magnetic field region. With $m_{0}=5$ MeV, we
obtain a flat curve for $T^{0}_{m}$ in the medium magnetic field region, and
in the weak and strong magnetic field region, $T^{0}_{m}$ decreases. With a
larger value of current quark mass, such as $m_{0}=10,\ 20,\ 50$ MeV,
$T^{0}_{m}$ is a monotonic decreasing functions of magnetic field. In the weak
magnetic field region, the Mott transition temperature is higher for larger
$m_{0}$, but in the strong magnetic field region, the Mott transition
temperature is lower for larger $m_{0}$. For $\pi^{+}$ meson, the Mott
transition temperature $T^{+}_{m}$ is only quantitatively modifies by current
quark mass effect. Associated with some oscillations, $T^{+}_{m}$ increases
fast in the weak magnetic field region, decreases in the medium magnetic field
region and slightly increases in the strong magnetic field region. With larger
$m_{0}$, $T^{+}_{m}$ becomes higher.

Due to the interaction with magnetic field, the charged pions $\pi^{\pm}$ show
different behavior from the neutral pion $\pi^{0}$. One common character of
pions $(\pi^{0},\ \pi^{\pm})$ is the mass jump at their Mott transitions,
which is induced by the dimension reduction of the constituent quarks, and
independent on the IMC effect and CQM effect. As a consequence of such jumps,
it may result in some interesting phenomena in relativistic heavy ion
collisions where the strong magnetic field can be created. For instance, when
the formed fireball cools down, there might be sudden enhancement/suppression
of pions.

Acknowledgement: Dr. Luyang Li is supported by Natural Science Basic Research
Plan in Shaanxi Province of China (Program No.2023-JC-YB-570). Prof. Shijun
Mao is supported by the NSFC Grant 12275204 and Fundamental Research Funds for
the Central Universities.

## References

  * (1) C. M. Ko, Z. G. Wu, L. H. Xia and G. E. Brown, Phys. Rev. Lett. 66, 2577 (1991) [Erratum-ibid. 67, 1811 (1991)].
  * (2) G. Q. Li and G. E. Brown, Phys. Rev. C 58, 1698 (1998)
  * (3) K. Paech, A. Dumitru, J. Schaffner-Bielich, H. Stoecker, G. Zeeb, D. Zschiesche and S. Schramm, Acta Phys. Hung. A 21, 151 (2004).
  * (4) J. Rafelski and B. Muller, Phys. Rev. Lett. 36, (1976) 517.
  * (5) D. E. Kharzeev, L. D. McLerran and H. J. Warringa, Nucl. Phys. A 803, (2008) 227.
  * (6) V. Skokov. A. Y. Illarionov and T. Toneev, Int. J. Mod. Phys. A 24, (2009) 5925.
  * (7) W. T. Deng and X. G. Huang, Phys. Rev. C 85, (2012) 044907; Phys. Lett. B 742, (2015) 296.
  * (8) K. Tuchin, Adv. High Energy Phys. (2013) 490495.
  * (9) N. Agasian and I. Shushpanov, J. High Energy Phys. 10 (2001) 006.
  * (10) G. Colucci, E. Fraga and A. Sedrakian, Phys. Lett. B 728, 19 (2014).
  * (11) J. Anderson, J. High Energy Phys. 10 (2012) 005, Phys. Rev. D 86, 025020 (2012).
  * (12) K. Kamikado and T. Kanazawa, J. High Energy Phys. 03 (2014) 009.
  * (13) G. Krein,and C. Miller, Symmetry, 13(4): 551 (2021).
  * (14) A. Ayala, J. L. Hern$\acute{a}$ndez, L. A. Hern$\acute{a}$ndez, R. L. S. Farias, R. Zamora, Phys. Rev. D 103, 054038 (2021).
  * (15) R. M. Aguirre, Phys. Rev. D 96, 096013 (2017).
  * (16) A. Ayala, R. L. S. Farias, S. Hernández-Ortiz, L. A. Hernández, D. M. Paret and R. Zamora, Phys. Rev. D 98, 114008 (2018).
  * (17) A. N. Tawfik, A. M. Diab and T. M. Hussein, Chin. Phys. C 43, 034103 (2019).
  * (18) A. Das and N. Haque, Phys. Rev. D 101, 074033 (2020).
  * (19) G. S. Bali, F. Bruckmann, G. Endrődi, Z. Fodor, S. D. Katz, S. Krieg, A. Schäfer and K. K. Szabó, J. High Energy Phys. 02 (2012) 044.
  * (20) Y. Hidaka and A. Yamatomo, Phys. Rev. D 87, 094502 (2013).
  * (21) E. Luschevskaya, O. Solovjeva, O. Kochetkov and O. Teryaev, Nucl. Phys. B 898, 627 (2015).
  * (22) G. S. Bali, B. Brandt, G. Endrődi and B. Gläßle, PoS LATTICE 2015, 265 (2016).
  * (23) E. Luschevskaya, O. Solovjeva and O. Teryaev, Phys. Lett. B 761, 393 (2016).
  * (24) G. S. Bali, B. Brandt, G. Endrődi and B. Gläßle, Phys. Rev. D 97, 034505 (2018).
  * (25) H. T. Ding, S. T. Li, A. Tomiya, X. D. Wang and Y. Zhang, Phys. Rev. D 104, 014505 (2021).
  * (26) H. T. Ding, S. T. Li, J. H. Liu, and X. D. Wang, Phys. Rev. D 105, 034514 (2022).
  * (27) S. Klevansky, Rev. Mod. Phys. 64, 649 (1992).
  * (28) S. Fayazbakhsh, S. Sadeghian and N. Sadooghi, Phys. Rev. D 86, 085042 (2012).
  * (29) S. Fayazbakhsh and N. Sadooghi, Phys. Rev. D 88, 065030 (2013).
  * (30) V. D. Orlovsky and Y. A. Simonov, J. High Energy Phys. 09 (2013) 136.
  * (31) R. Zhang, W. J. Fu and Y. X. Liu, J. Eur. Phys. C 76, 307 (2016).
  * (32) S. Avancini, W. Travres and M. Pinto, Phys. Rev. D 93, 014010 (2016).
  * (33) K. Hattori, T. Kojo and N. Su, Nucl. Phys. A 951, 1 (2016).
  * (34) Y. A. Simonov, Phys. Atom. Nucl. 79, 455 (2016). [Yad. Fiz. 79, 277 (2016)].
  * (35) M. A. Andreichikov, B. O. Kerbikov, E. V. Luschevskaya, Y. A. Simonov and O. E. Solovjeva, J. High Energy Phys. 05 (2017) 007.
  * (36) S. Avancini, R. Farias, M. Pinto, W. Travres and V. Tim$\acute{o}$teo, Phys. Lett. B 767, 247 (2017).
  * (37) S. J. Mao and Y. X. Wang, Phys. Rev. D 96, 034004 (2017).
  * (38) M. A. Andreichikov and Y. A. Simonov, Eur. Phys. J. C 78, 902 (2018).
  * (39) C. A. Dominguez, M. Loewe and C. Villavicencio, Phys. Rev. D 98, 034015 (2018).
  * (40) Z. Y. Wang and P. F. Zhuang, Phys. Rev. D 97, 034026 (2018).
  * (41) M. Coppola, D. Dumm and N. Scoccola, Phys. Lett B 782, 155 (2018).
  * (42) H. Liu, X. Wang, L. Yu and M. Huang, Phys. Rev. D 97, 076008 (2018).
  * (43) D. G. Dumm, M. I. Villafañe and N. N. Scoccola, Phys. Rev. D 97, 034025 (2018).
  * (44) S. S. Avancini, R. L. S. Farias and W. R. Tavares, Phys. Rev. D 99, 056009 (2019).
  * (45) N. Chaudhuri, S. Ghosh, S. Sarkar and P. Roy, Phys. Rev. D 99, 116025 (2019).
  * (46) M. Coppola, D. G. Dumm, S. Noguera and N. N. Scoccola, Phys. Rev. D 100, 054014 (2019).
  * (47) J. Y. Chao, Y. X. Liu and L. Chang, arXiv: 2007.14258.
  * (48) S. J. Mao, Phys. Rev. D 99, 056005 (2019).
  * (49) J. Moreira, P. Costa, and T. E. Restrepo, Phys. Rev. D 102, 014032 (2020).
  * (50) D. N. Li, G. Q. Cao and L. Y. He, Phys. Rev. D 104, 074026 (2021).
  * (51) B. K. Sheng, Y. Y. Wang, X. Y. Wang, L. Yu, Phys. Rev. D 103, 094001 (2021).
  * (52) S. J. Mao, Chin. Phys. C 45, 021004 (2021).
  * (53) K. Xu, J.Y. Chao and M. Huang, Phys. Rev. D 103, 076015 (2021).
  * (54) S. S. Avancini, M. Coppola, N. N. Scoccola, and J. C. Sodré, Phys. Rev. D 104, 094040 (2021).
  * (55) S. J. Mao and Y. M. Tian, Phys. Rev. D 106, 094017 (2022).
  * (56) L. Y. Li and S. J. Mao, Chin. Phys. C 46, 094105 (2022).
  * (57) J. P. Carlomagno, D. G. Dumm, M. F. I. Villafa$\tilde{n}$e, S. Noguera, and N. N. Scoccola, Phys. Rev. D 106, 094035 (2022).
  * (58) J. P. Carlomagno, D. G. Dumm, S. Noguera, and N. N. Scoccola, Phys. Rev. D 106, 074002 (2022).
  * (59) J. Mei and S.J. Mao, Phys. Rev. D 107, 074018 (2023).
  * (60) G. S. Bali, F. Bruckmann, G. Endr$\ddot{o}$di, Z. Fodor, S. D. Katz, S. Krieg, A. Sch$\ddot{a}$fer and K. K. Szab$\acute{o}$, Phys. Rev. D 86, 071502 (2012).
  * (61) G. S. Bali, F. Bruckmann, G. Endr$\ddot{o}$di, F. Gruber and A. Sch$\ddot{a}$fer, J. High Energy Phys. 04, 130(2013).
  * (62) V. Bornyakov, P.V. Buividovich, N. Cundy, O. A. Kochetkov and A. Sch$\ddot{a}$fer, Phys. Rev. D 90, 034501 (2014).
  * (63) E. M. Ilgenfritz, M. Muller-Preussker, B. Petersson and A. Schreiber, Phys. Rev. D 89, 054512 (2014).
  * (64) G. Endr$\ddot{o}$di, J. High Energy Phys. 07, 173 (2015).
  * (65) M. D’Elia, F. Manigrasso, F. Negro and F. Sanfilippo, Phys. Rev. D 98, 054509 (2018).
  * (66) G. Endr$\ddot{o}$di, M. Giordano, S. D. Katz, T. G. Kovacs and F. Pittler, J. High Energy Phys. 07, 009 (2019).
  * (67) A. Tomiya, H.T. Ding, X.D. Wang, Y. Zhang, S. Mukherjee and C. Schmidt, PoS Lattice2018, 163 (2019).
  * (68) E. S. Fraga and A. J. Mizher, Phys. Rev. D 78, 025016 (2008).
  * (69) C. V. Johnson and A. Kundu, J. High Energy Phys. 12, 053 (2008).
  * (70) A. J. Mizher, M. N. Chernodub and E. S. Fraga, Phys. Rev. D 82, 105016 (2010).
  * (71) K. Fukushima, M. Ruggieri and R. Gatto, Phys. Rev. D 81, 114031 (2010).
  * (72) R. Gatto and M. Ruggieri, Phys. Rev. D 82, 054027 (2010); Phys. Rev. D 83, 034016 (2011).
  * (73) F. Preis, A. Rebhan and A. Schmitt, J. High Energy Phys. 03, 033 (2011).
  * (74) V. Skokov, Phys. Rev. D 85, 034026 (2012).
  * (75) E. S. Fraga, J. Noronha and L. F. Palhares, Phys. Rev. D 87, 114014 (2013).
  * (76) K. Fukushima and Y. Hidaka, Phys. Rev. Lett 110, 031601 (2013).
  * (77) J. Y. Chao, P. C. Chu and M. Huang, Phys. Rev. D 88, 054009 (2013).
  * (78) F. Bruckmann, G. Endr$\ddot{o}$di and T. G. Kovacs, arXiv: 1311.3178.
  * (79) T. Kojo and N. Su, Phys. Lett. B 720, 192 (2013).
  * (80) F. Bruckmann, G. Endrodi and T. G. Kovacs, J. High Energy Phys. 04, 112 (2013).
  * (81) A. Ayala, M. Loewe, A. J$\acute{u}$lia Mizher and R. Zamora, Phys. Rev. D 90, 036001 (2014).
  * (82) A. Ayala, L. Alberto Hern$\acute{a}$ndez, A. J$\acute{u}$lia Mizher, J. Crist$\acute{o}$bal Rojas and C. Villavicencio, Phys. Rev. D 89, 116017 (2014).
  * (83) R. L. S. Farias, K. P. Gomes, G. Krein and M. B. Pinto, Phys. Rev. C 90, 025203 (2014).
  * (84) M. Ferreira, P. Costa, O. Lourenco, T. Frederico and C. Provid$\hat{e}$ncia, Phys. Rev. D 89, 116011 (2014).
  * (85) M. Ferreira, P. Costa and C. Provid$\hat{e}$ncia, Phys. Rev. D 89, 036006 (2014).
  * (86) M. Ferreira, P. Costa, D. P. Menezes, C. Provid$\hat{e}$ncia and N. N. Scoccola, Phys. Rev. D 89, 016002 (2014).
  * (87) P. Costa, M. Ferreira, H. Hansen, D. P. Menezes and C. Provid$\hat{e}$ncia, Phys. Rev. D 89, 056013 (2014).
  * (88) E. J. Ferrer, V.de la Incera, I. Portillo and M. Quiroz, Phys. Rev. D 89, 085034 (2014).
  * (89) A. Ayala, C. A. Dominguez, L. A. Hern$\acute{a}$ndez, M. Loewe and R. Zamora, Phys. Rev. D 92, 096011 (2015).
  * (90) N. Mueller and J.M. Pawlowski, Phys. Rev. D 91, 116010 (2015).
  * (91) E. J. Ferrer, V.de la Incera and X. J. Wen, Phys. Rev. D 91, 054006 (2015).
  * (92) J. Braun, W. A. Mian and S. Rechenberger, Phys. Lett. B 755, 265 (2016).
  * (93) S. J. Mao, Phys. Lett. B 758, 195 (2016); Phys. Rev. D 94, 036007 (2016); Phys. Rev. D 97, 011501(R) (2018).
  * (94) J. Mei and S.J. Mao, Phys. Rev. D 102, 114035 (2020).
  * (95) A. Ayala, L. A. Hern$\acute{a}$ndez, M. Loewe and C. Villavicencio, Eur. Phys. J. A 57, 234 (2021).
  * (96) S. J. Mao, Phys. Rev. D 106, 034018 (2022).
  * (97) Y. Nambu and G. Jona-Lasinio, Phys. Rev. 122, 345 (1961) and 124, 246 (1961).
  * (98) M. K. Volkov, Phys. Part. Nucl. 24, 35 (1993).
  * (99) T. Hatsuda and T. Kunihiro, Phys. Rep. 247, 221 (1994).
  * (100) M. Buballa, Phys. Rep. 407, 205 (2005).
  * (101) P. Zhuang, J. H$\ddot{u}$fner, S. P. Klevansky, Nucl. Phys. A 576, 525 (1994).
  * (102) J. Goldstone, Nuovo Cim. 19, 154-164 (1961).
  * (103) J. Goldstone, A. Salam and S. Weinberg, Phys. Rev. 127, 965-970 (1962).
  * (104) N. F. Mott, Rev. Mod. Phys. 40, 677 (1968).
  * (105) J. Huefner, S. Klevansky, and P. Rehberg, Nucl. Phys. A 606, 260 (1996).
  * (106) P. Costa, M. Ruivo, and Y. Kalinovsky, Phys. Lett. B 560, 171 (2003).
  * (107) H. Liu, L. Yu, M. Chernodub and M. Huang, Phys. Rev. D 94, 113006 (2016).
  * (108) A. Ayala, C. A. Dominguez, L. A. Hern$\acute{a}$ndez, M. Loewe, A. Raya, J. C. Rojas and C. Villavicenico, Phys. Rev. D 94, 054019 (2016).
  * (109) R. L. S. Farias, V. S. Timoteo, S. S. Avancini, M. B. Pinto and G. Klein, Eur. Phys. J. A 53, 101 (2017).