\section{Introduction}
\setcounter{footnote}{0}
\label{intro}
The North American Nanohertz Observatory for Gravitational Waves ({\rm NANOGrav}~) searches for an isotropic stochastic gravitational wave background (SGWB) by analyzing the cross-power spectrum of pulsar timing residuals~\cite{Arzoumanian:2018saf}. The recent $12.5$ year pulsar timing array (PTA) data released by the {\rm NANOGrav}~ has reported the discovery of a stochastic common-spectrum process~\cite{Arzoumanian:2020vkk}, which can be fitted into a power law in a narrow range of frequencies. At this moment, they report that the angular correlations are inconclusive and therefore, the collaboration does not claim a positive signal of gravitational waves (GWs). However, if the signal is interpreted as an upper bound on GW background, we can explore the primordial universe with the {\rm NANOGrav}~ constraints.
Assuming the {\rm NANOGrav}~ signal to be a GW background, the source of such a signal can be supermassive black hole merger events~\cite{Sesana:2004sp,Arzoumanian:2020vkk}, cosmic strings~\cite{Ellis:2020ena,Blasi:2020mfx,Buchmuller:2020lbh,Samanta:2020cdk}, phase transition of a hidden sector\cite{Nakai:2020oit,Addazi:2020zcj,Ratzinger:2020koh,Li:2020cjj}, inflationary perturbations in the first~\cite{Vagnozzi:2020gtf} or second orders~\cite{Vaskonen:2020lbd,DeLuca:2020agl,Kohri:2020qqd,Domenech:2020ers} etc.~\cite{Pandey:2020gjy}.

The LIGO/VIRGO detection of GWs from the binary black hole mergers has opened up a new window towards the early universe cosmology~\cite{Abbott:2016blz,Abbott:2016nmj,TheLIGOScientific:2016pea,Abbott:2016bqf,Abbott:2017vtc,Abbott:2017oio,Abbott:2017gyy}. The probability that the observed black holes of supersolar masses are of primordial origin has also rekindled the interest in the primordial black holes (PBHs)~\cite{Fernandez:2019kyb}. Other than such merger events, GWs of a stochastic nature can also originate from the primordial tensor perturbations generated during inflation~\cite{Caprini:2018mtu,Christensen:2018iqi}. The amplitude of the primordial tensor power spectrum at the cosmic microwave background (CMB) pivot scale ($k_{\rm pivot}=0.05$ Mpc$^{-1}$) is constrained by the Planck 2018 and BICEP2/Keck Array~\cite{Ade:2018gkx,Aghanim:2019ame} observations in terms of the tensor-to-scalar ratio $r<0.056$ with $95\%$ confidence limit. \pri{However, if some feature in the inflationary dynamics induces a blue-tilt in the tensor power spectrum at smaller scales, GWs of larger amplitude are possible~\cite{Cook:2011hg,Barnaby:2011qe,Wang:2014kqa,Kuroyanagi:2014nba}. Since the {\rm NANOGrav}~ observation are at much higher frequencies, $f\in (2.5\times10^{-9},1.2\times10^{-8})$, as compared to the CMB frequencies, a blue tilted tensor power spectrum may be consistent with both the CMB and {\rm NANOGrav}~ observations.}

However, in the second and higher orders of perturbation theory, GWs are sourced by the first order curvature perturbations~\cite{Ananda:2006af,Baumann:2007zm,Espinosa:2018eve,Kohri:2018awv}.
If the curvature perturbations have a red-tilted power spectrum throughout inflation, such second order GWs are suppressed with respect to the first order contribution. \pri{Observation of the temperature and polarization anisotropies in the CMB sky has constrained the primordial scalar power spectrum to have a red tilt ($n_s=0.9649$) and tiny amplitude $A_s \sim \mathcal{O}(10^{-9})$ at the CMB pivot scale~\cite{Ade:2018gkx,Aghanim:2019ame}. However, if there is a large running of $n_s$, the primordial scalar power spectrum can have a blue-tilt at smaller scales, and then the primordial curvature perturbations can be large at these scales, which can induce second order GW of large amplitude (see~\cite{Clesse:2015wea,Garcia-Bellido:2017mdw,Germani:2017bcs,Ballesteros:2017fsr,Bhaumik:2019tvl,Mishra:2019pzq,Ballesteros:2020qam,Palma:2020ejf,Fumagalli:2020adf,Braglia:2020eai,Ozsoy:2020kat,Aldabergenov:2020yok,Braglia:2020taf} for relevant physical models). Such large curvature perturbations, leading to density contrast larger than the critical value for gravitational collapse, will lead to a copious amount of PBH production.}
 Hence, the phenomenologies of PBH and GW in the second order are always related~\cite{Hawking:1971ei,Carr:1974nx,Carr:1975qj}. PBH abundance with the {\rm NANOGrav}~ result has been studied in~\cite{Vaskonen:2020lbd,DeLuca:2020agl,Kohri:2020qqd,Domenech:2020ers}.

Big Bang Nucleosynthesis (BBN) constraints demand that the universe was radiation dominated at least at the temperature $T_{\rm BBN}\simeq 5$ MeV. CMB constraint on $r$ reveals that the upper bound on the energy scale of inflaton potential is $V\simeq 10^{16}$ GeV~\cite{Akrami:2018odb,Akrami:2019izv}.
For a standard slow roll-inflation, the energy scale of inflation remains nearly constant until the end of inflation. The epoch between the end of inflation and BBN, which can span a large number of e-folds, is not well constrained by the observations. In the standard theoretical description, the universe quickly becomes radiation dominated after the end of inflation, and the standard radiation epoch sustains until the matter radiation equality ($T_{\rm eq}$). However, it is possible that the equation of state $w$ of the universe deviates from radiation after the end of inflation and before BBN\footnote{For a review, see~\cite{Allahverdi:2020bys}.}.
An epoch of nonstandard post inflationary evolution can be motivated by slow reheating ($w$ slowly changes from $0$ to $1/3$)~\cite{Carr:2018nkm}; early matter dominated epoch ($w=0$) with energy dominated by the moduli field~\cite{Coughlan:1983ci,Banks:1993en,deCarlos:1993wie,Conlon:2005jm,Cicoli:2016olq,Bhattacharya:2017ysa,Maharana:2017fui}; kinetic energy domination ($w=1$) in case of the quintessential inflation models~\cite{Peebles:1998qn}; general stiff domination~\cite{DiMarco:2018bnw} etc. The dynamics of formation and evolution of PBH in such epochs is also different and their abundance is related to the amplitude of the second order GW. Therefore, any deviation from the standard radiation epoch before BBN affects the PBH and GW dynamics.

In this paper, we inspect the first and second order primordial GWs in a non-standard post-inflationary epoch (i.e., equation of state $w\neq 1/3$) that is active at the scales probed by the {\rm NANOGrav}~ . \suk{The evolution of the first order tensor perturbations upon horizon entry depends on the background as $\propto a^{3(1+w)-4}$. Thus, for a standard radiation dominated (RD) epoch, the scale dependence of the first order GW spectrum is only due to the primordial tensor tilt $n_t$. However, for $w\neq 1/3$, the slope of the resulting GW spectra depends on $w$. Here, we study the first order GW spectra for $w=1/3$, $0$, and $1$ and find that only a fraction of the parameter space for the $w=0$ scenario is consistent with the {\rm NANOGrav}~ data, when the constraints from CMB, BBN, and LIGO are included. We compare our analysis for the first order GW spectrum in RD with ref.~\cite{Vagnozzi:2020gtf}.

The critical value $\delta _c$ of the density contrast strongly depends on the background in which the scalar perturbations enter the horizon and collapse to form PBHs. This has motivated many studies of PBH formation in nonstandard epochs~\cite{1981SvA....25..406P,1982AZh....59..639P,Hwang:2012bi,Alabidi:2013lya,Harada:2016mhb,Carr:2017edp,Harada:2017fjm,Inomata:2019zqy,Inomata:2019ivs,Figueroa:2019paj,Matsubara:2019qzv,Bhattacharya:2019bvk}. In this paper, to analyse second order GW spectra for abundant PBH formation, we have considered two nonstandard epochs: $w=1$, i.e. a kination epoch and $w=1/9$, which is a dustlike epoch. Since $\delta _c$ for both of these nonstandard epochs are smaller than $\delta _c^{\rm RD}$ (see Fig.3 in ref.~\cite{Harada:2013epa}), we expect to reach a given PBH abundance with much lower amplitudes of the primordial scalar perturbations than that required in RD. This work focuses on determining whether those values of the primordial scalar amplitude can explain the {\rm NANOGrav}~ observation. We follow the formation of PBH and evolution of the scalar induced GW spectra in second order in $w=1$ and $w=1/9$ epochs for a few functional forms of the primordial curvature power spectrum. We find that although abundant PBH formation in a kination epoch is not consistent with the {\rm NANOGrav}~ detection, for the dustlike epoch, large PBH abundance can explain the {\rm NANOGrav}~ data as a second order GW spectrum. The standard case of PBH formation and second order GW evolution in RD has also been studied here and compared with the contemporary works with a similar approach. We find our results for the RD epoch to be consistent with ref.~\cite{Kohri:2020qqd}, however, disagree with those of ref.~\cite{Vaskonen:2020lbd} due to reasons discussed later.}

The rest of the paper is organised as follows: in Sec.~\ref{GWsec}, we represent the result from the {\rm NANOGrav}~ in terms of the theoretical quantity: the GW spectrum $\Omega_{\rm GW}$. In Sec.~\ref{GW1sec}, we discuss the effect of a general equation of state $0\leq w\leq 1$ on the first order GW and constrain the relevant parameters using the power law nature of the observed GW signal, fitted by {\rm NANOGrav}~ using the frequency in the first five  bins. Sec.~\ref{GW2sec} relates the second order GW to the observation and translates the constraints in terms of PBH abundance. In Sec.~\ref{ressec}, we discuss our results and conclude.

\section{NANOGrav Result and Gravitational waves}
\label{GWsec}
PTA experiments typically describe the obtained result for GW background in terms of the characteristic strain spectrum $h_c(f)$. The quantity $h_c(f)$ is typically fitted with a power law dependence on frequency $f$~\cite{Zhao:2013bba,Liu:2015psa,Arzoumanian:2018saf}
\begin{eqnarray}
h_c(f) = A_{\rm GWB} \left ( \frac{f}{f_{\rm yr}} \right )^{\alpha_{_{\rm GWB}}}\,,
\label{eq:strain}
\end{eqnarray}
with $\alpha_{_{\rm GWB}}=(3-\gamma)/2$, where $\gamma$ is the timing-residual cross-power spectral density. The reference frequency is $f_{\rm yr} = {\rm yr}^{-1}=3.1\times 10^{-8}$ Hz.
The recent 12.5 yr observation from the {\rm NANOGrav}~~\cite{Arzoumanian:2020vkk} measured
 the characteristic strain in the frequency range $f\in (2.5\times10^{-9},1.2\times10^{-8})$ Hz and found the evidence of a stochastic common-spectrum process (CP). The data is fitted as a power law and the corresponding GW energy density $\Omega_{\rm GW}(f)$ is given as~\cite{Zhao:2013bba,Liu:2015psa,Arzoumanian:2018saf}:
\begin{eqnarray}
\Omega_{\rm GW}(f) = \frac{2\pi^2}{3H_0^2}f^2h_c^2(f)=\frac{2\pi^2f_{\rm yr}^2}{3H_0^2}A_{\rm CP}^2\bigg(\frac{f}{f_{\rm yr}}\bigg)^{5-\gamma _{\rm CP}}\,.
\label{eq:hcGW}
\end{eqnarray}
PTA data from {\rm NANOGrav}~ give a joint $A_{_{\rm CP}}$-$\gamma_{_{\rm CP}}$ posterior distribution.
Therefore, early universe theories with a prediction for a stochastic GW signal can be constrained using this result in terms of the amplitude $A_{_{\rm CP}}$ and slope $\gamma_{_{\rm CP}}$.

In a perturbed spatially flat Friedmann-Lema\^{i}tre-Robertson-Walker (FLRW) universe, the tensor perturbations $h_{ij}$ are traceless and transverse. The tensor modes are given as:
\begin{equation}
h_{ij} (\eta,{\bf x}) = \sum_{\lambda} \int  \frac{d^3 {\bf k}}{(2\pi)^{3}} e^{i {\bf k.x}} \epsilon_{ij}^\lambda({\bf k})h_{{\bf k}}^\lambda({\eta}) ,
\end{equation}
where $\eta$ is conformal time, ${\bf k}$ is the comoving wavenumber and $ \epsilon_{ij}^\lambda({\bf k})$ are the polarization tensors for the two polarization states $\lambda = + ,\times $. The tensor power spectrum $\mathcal{P}_h$ is defined as:
\begin{equation}
\frac{k^3}{2\pi^{2}} \langle h_{{\bf k}}^\lambda({\eta})h_{{\bf k'}}^{\lambda'}({\eta}) \rangle = \delta_{\lambda \lambda'} \delta^3({\bf k}+{\bf k'})\mathcal{P}_h(k,\eta)\,.
\end{equation}
The GW energy density spectrum $\Omega_{{\rm GW}}$ is defined as the GW energy density $\rho_{{\rm GW}}$ in a comoving wavenumber interval $(k, k + dk)$, normalised with the critical density $\rho_{c}(\eta)$,
\begin{equation}
\label{originalGWPh}
\Omega_{{\rm GW}} (\eta ,k) \equiv \frac{\rho_{\rm GW}(\eta,k)}{\rho_c (\eta)} = \frac{1}{24}\bigg(\frac{k}{ \mathcal{H}}\bigg)^2 \mathcal{P}_h(k,\eta)\, .
\end{equation}
Observationally relevant quantity is the GW spectrum at present $\Omega_{{\rm GW}}(\eta_0 ,k)$.
As discussed in Sec.~\ref{intro}, for a given theory, GW spectrum can be evaluated at first and higher orders of perturbation theory. Thus, theoretical predictions for GW resulting in different orders can be constrained with the {\rm NANOGrav}~ observation. When the scalar fluctuations give rise to PBH, then the PBH abundance can also be probed in terms of the observed GW spectra, \suk{where the second order contribution is dominant over the first order.}

In this work, we focus on analysing the first and second order GW in a nonstandard cosmological history, where the background evolution deviates from the standard radiation domination (RD) before BBN. We separately derive the predictions for different orders and compare them with the {\rm NANOGrav}~ observation to discuss bounds on the relevant primordial quantities.

\section{Primordial gravitational waves in first order}
\label{GW1sec}
Inflation predicts a stochastic GW background arising from the tensor fluctuations~\cite{Caprini:2018mtu,Christensen:2018iqi}. \pri{In the absence of anisotropic stresses, the evolution of these tensor fluctuations is source free at the first order of perturbations.} GW background depends on the primordial tensor power spectrum, which can be written as:
\begin{equation}
\mathcal{P}_h(k,\eta)= T(k,\eta) \Delta ^2_{h,{\rm inf}}(k),\label{omtot}
\end{equation}
where $\Delta ^2_{h,{\rm inf}}(k) =\frac{2}{\pi^{2}}\bigg(\frac{H_{\rm inf}}{M_{\rm P}}\bigg)^2 \bigg(\frac{k}{k_{p}}\bigg)^{n_t}$ is the primordial tensor power spectrum at horizon re-entry and $T(k,\eta)$ is the transfer function. Here $n_t$ is the spectral index, $k_{p}$ is a pivot scale and $H_{\rm inf}$ is the Hubble rate when mode $k_{p}$ exited the horizon during inflation. For a mode $k$ re-entering the horizon at time $\eta(k)$, the transfer function depends on the background evolution from the time of its horizon re-entry until the time of its observation.
In the case of standard cosmological evolution history, radiation domination (RD) starts just after inflation ends (assuming instant reheating) followed by the matter domination (MD). However, there can be a pre-BBN epoch with equation of state (EoS) other than the radiation ($w\neq1/3$) in alternative cosmological histories. In such alternative cosmological histories, the transfer function will be different for the modes re-entering the horizon during the nonstandard epoch. In this work, we have assumed a period with EoS $w\neq1/3$ from the end of inflation till the onset of RD at temperature $T_1$. GW energy density in such a nonstandard cosmological thermal history at the first order of perturbations is \cite{Bernal:2019lpc,Figueroa:2019paj,Bernal:2020ywq}:
\begin{eqnarray}
\Omega_{\rm GW,0}^{(1)}(k)&=&\frac{\Omega_{\rm rad, 0}}{12\pi^2}\bigg(\frac{g_{*,k}}{g_{s,k}}\bigg)\bigg(\frac{g_{s,0}}{g_{s,k}}\bigg)^{4/3}\bigg(\frac{H_{\rm inf}}{M_{\rm P}}\bigg)^2\frac{\Gamma ^2 (\alpha +1/2)}{2^{2(1-\alpha}\alpha ^{2\alpha}\Gamma ^2(3/2)}\mathcal{W}(\kappa)\kappa^{2(1-\alpha)} \bigg(\frac{k}{k_{p}}\bigg)^{n_t},\label{GW_1_gen}
\end{eqnarray}
where $\alpha=\frac{2}{1+3w}$, $\kappa=\frac{k}{k(T_1)}=\frac{f}{f(T_1)}$, and
\begin{equation}
\mathcal{W}(\kappa)= \frac{\pi \alpha}{2\kappa}\bigg[\bigg(\kappa J_{\alpha+1/2}(\kappa)-J_{\alpha-1/2}(\kappa)\bigg)^2+\kappa^2J^2 _{\alpha-1/2}(\kappa)\bigg],\label{Bessel_GW1}
\end{equation}
\pri{where $J_{i}$ and $\Gamma(\alpha)$ are the Bessel function of order $i$ and the gamma function respectively. $k(T_1)$ and $f(T_1)$ are the wavenumber and frequency corresponding to the temperature ($T_1$) of the universe at the time of the onset of RD epoch respectively. } Also, $\Omega_{\rm rad, 0}=9\times10^{-5}$ is the radiation energy fraction in the present universe, $g_{*}$ and $g_{s}$ are the number of relativistic degrees of freedom for the energy and entropy respectively. We also find that $\mathcal{W}(\kappa) \approx \alpha$ for $f> f(T_1)$. It is clear from Eq.\eqref{GW_1_gen} that the GW spectrum today depends on $H_{\rm inf}$, $w$, $n_t$, and $f(T_1)$ (for detailed discussion see~ref~\cite{Bernal:2019lpc,Figueroa:2019paj,Bhattacharya:2019bvk}). First order GW spectrum for the three different values of EoS,  $w=0,1/3$, and $1$, are shown in Fig.~\ref{fig:first_GW_kd_rd_md} where we have fixed $n_t = 0$, $H_{\rm inf} = 6.2 \times 10^{13}$ GeV and $f(T_1)=10^{-10}$ Hz. In order to do the analysis with the {\rm NANOGrav}~ result, it is useful to relate the {\rm NANOGrav}~ observational parameters $A_{_{\rm CP}}$ and $\gamma_{_{\rm CP}}$ to the set of theoretical parameters: ($H_{\rm inf}$, $w$, $n_t$, $f(T_1)$). Using the Eq.s~\eqref{eq:strain},  \eqref{eq:hcGW} and \eqref{GW_1_gen}, we get
\begin{eqnarray}
 A_{\rm CP} = \bigg( \frac{3 H_0^2}{2 \pi^2}\frac{\Omega _{\rm rad, 0}}{12\pi ^2}\bigg(\frac{g_{*,k}}{g_{s,k}}\bigg)\bigg(\frac{g_{s,0}}{g_{s,k}}\bigg)^{4/3} \frac{\Gamma ^2 (\alpha +1/2)}{2^{2(1-\alpha}\alpha ^{2\alpha}\Gamma ^2(3/2)} \alpha  \bigg)^{1/2}\frac{H_{\rm inf}}{M_{\rm P}} \bigg(
 \frac{f_{\rm yr}}{f_{p}}\bigg)^{n_t/2} \frac{{f_{\rm yr}}^{\alpha}}{{f_{T1}}^{1-\alpha}} \,, \label{eq:acp}
\end{eqnarray}
and
\begin{eqnarray}
\gamma_{_{\rm CP}} = 3+2\alpha-n_t\,.
\label{eq:gammacp}
\end{eqnarray}
\begin{figure}[!b]
\begin{center}
\includegraphics[width=0.6\textwidth]{first_order_GW_kd_rd_md_with_obs.pdf}
\caption{ \small{First order GW spectrum for the three different cosmological histories are shown here. We have fixed $n_t = 0$, $H_{\rm inf} = 6.2 \times 10^{13}$ GeV and $f(T_1)=10^{-10}$ Hz. LIGO runs for O1,O2 and O5 (future) are given as cyan curves (top to bottom); magenta shaded region represnents the proposed sensitivity for the LISA survey; black shaded region represents the $1\sigma$ bound from the {\rm NANOGrav}~ 12.5 year observation; orange shaded part represents the proposed sensitivity from the future SKA survey. Purple dashed line shows the GW spectrum in cosmologcal history with an early MD epoch and $n_t=1$. }}\label{fig:first_GW_kd_rd_md}
\end{center}
\end{figure}
Evidently, Eq.s~\eqref{eq:acp} and \eqref{eq:gammacp} show that $A_{\rm CP}$ depends on $H_{\rm inf}$, $w$, $n_t$, and $f(T_1)$, whereas $\gamma _{\rm CP}$ depends only on $w$ and $n_t$.
\begin{figure}[!t]
  \subfloat[\label{fig:acp_w}]{
	\begin{minipage}[c][1\width]{
	   0.5\textwidth}
	   \centering
	   \includegraphics[width=1.05\textwidth]{Acp_vs_w_fixed_hinf_nt_ft1.pdf}
	\end{minipage}}
 \hfill
  \subfloat[\label{fig:gammacp_w}]{
	\begin{minipage}[c][1\width]{
	   0.5\textwidth}
	   \vspace{0.3cm}
	   \centering
	   \includegraphics[width=0.98\textwidth]{gammacp_vs_w_fixed_nt_ft1.pdf}
	\end{minipage}}
\caption{\small{(a) $A_{{\rm CP}}$ is plotted against EoS $w$ for $H_{\rm inf} =10^{11}$ and $10^{13}$ GeV, $n_t=1,0,-0.1$. (b) $\gamma _{{\rm CP}}$ is plotted against EoS $w$ for $n_t=1,0,-0.1$. Gray shaded band in both the plots is the allowed region from the NANOGrav observation. }}
\end{figure}
In Fig.~\ref{fig:acp_w}, $A_{\rm CP}$ is plotted against EoS $w$ for $H_{\rm inf} =10^{11}$ and $10^{13}$ GeV, $n_t=1,0,-0.1$, and $f(T_1)=10^{-10}$ Hz, which corresponds to temperature $T_1=10$ MeV. We have also plotted $\gamma _{\rm CP}$ against EoS $w$ for $n_t=1,0,-0.1$ in Fig.~\ref{fig:gammacp_w}. Gray shaded bands in Fig.s~\ref{fig:acp_w} and~\ref{fig:gammacp_w} represent the $1\sigma$ and $2\sigma$ region allowed from the NANOGrav 12.5-yr result.
GWs generated during inflation leave their imprint on the CMB temperature  and polarisation anisotropy, which put a constraint on tensor-to-scalar ratio $r<0.056$ at $95\%$ confidence limit at the CMB pivot scale $k_{\rm pivot}=0.05$ Mpc$^{-1}$ from the joint analysis of BICEP2/Keck Array and Planck CMB data\cite{Akrami:2018odb}. This bound, in turn, puts a constraint on the value of $H_{{\rm inf}} < 6.2 \times 10^{13}$ GeV at $95\%$ confidence limit~~\cite{Akrami:2018odb,Akrami:2018vks,Aghanim:2018eyx,Aghanim:2019ame,Ade:2018gkx}.
Therefore, it is evident from  Fig.~\ref{fig:acp_w} that a positive value of $n_t$ is required to explain the allowed $A_{\rm CP}$ values. On the other hand, it can be seen from Fig.~\ref{fig:gammacp_w} that $\gamma_{\rm CP}$  value decreases with increasing $n_t$ and $w$ values. Therefore, for large values of $w$, $\gamma_{\rm CP}$ goes outside of the allowed region from the {\rm NANOGrav}~ observation.
\begin{figure}[]
  \subfloat[\label{fig:Hinf_nt_allowed_rd_mdw1}]{
	\begin{minipage}[c][1\width]{
	   0.5\textwidth}
	   \centering
	   \includegraphics[width=0.95\linewidth,height=0.8\linewidth]{allowed_hinf_nt_w1.pdf}
	\end{minipage}}
 \hfill
  \subfloat[\label{fig:Hinf_nt_allowed_rd_mdw13}]{
	\begin{minipage}[c][1\width]{
	   0.5\textwidth}
	   \centering
	   \includegraphics[width=0.95\linewidth,height=0.8\linewidth]{allowed_hinf_nt_w13.pdf}
	\end{minipage}}\\
	\subfloat[\label{fig:Hinf_nt_allowed_rd_mdw0}]{
	\begin{minipage}[c][1\width]{
	   0.5\textwidth}
	   \centering
	   \includegraphics[width=0.95\linewidth,height=0.8\linewidth]{allowed_hinf_nt_w0.pdf}
	\end{minipage}}
 \hfill
  \subfloat[\label{fig:Hinf_nt_allowed_rd_mdw0acp}]{
	\begin{minipage}[c][1\width]{
	   0.5\textwidth}
	   \centering
	   \includegraphics[width=0.95\linewidth,height=0.8\linewidth]{allowed_NanoGrav_param_w0.pdf}
	\end{minipage}}
\caption{ \small{ \pri{Allowed $H_{{\rm inf}}-n_t$ parameter space required to explain the NANOGrav result by the first order GW for different cases of cosmological thermal histories along with the bounds from different observations are shown here. Magenta shaded region represent the $H_{\rm inf}$ values excluded by the Planck CMB and BICEP2 observations. Green and yellow shaded regions represent the parameter space excluded by the LIGO constraints and by the contribution of GWs to $\Delta N_{\rm eff}$ respectively. (a) Gray shaded contours show the 1-$\sigma $ and 2-$\sigma$ region of $H_{{\rm inf}}-n_t$ parameter space required to explain the NANOGrav result by the first order GW in a pre-BBN kination dominated era. (b) Orange shaded contours represent the 1-$\sigma $ and 2-$\sigma$ region of $H_{{\rm inf}}-n_t$ parameter space required to explain the NANOGrav result by the first order GW in the standard cosmological thermal history. (c)  Blue shaded contours show the 1-$\sigma $ and 2-$\sigma$ region of $H_{{\rm inf}}-n_t$ parameter space required to explain the NANOGrav result by the first order GW in a pre-BBN matter dominated era.
(d) 1-$\sigma $ and 2-$\sigma$ contours of $\gamma_{{\rm CP}}-A_{{\rm CP}}$ values as reported by the {\rm NANOGrav}~ obvservation are shown here in black dashed and solid lines respectively. Blue shaded region shows the $\gamma_{{\rm CP}}-A_{{\rm CP}}$ parameter space which corresponds to the {\rm NANOGrav}~ allowed $H_{{\rm inf}}-n_t$ parameter space and consistent with the CMB, BBN and LIGO bounds for the nonstandard cosmological history with a pre-BBN MD era.}}}\label{fig:Hinf_nt_allowed_rd_md}
\end{figure}

\pri{In this work, our aim is to find the allowed region of the $H_{{\rm inf}}-n_t$ parameter space required to explain the NANOGrav observation for different values of the EoS. It is clear from Eq.s~\eqref{eq:acp} and \eqref{eq:gammacp} that there are four parameters $H_{\rm inf}$, $w$, $n_t$, and $f(T_1)$ which affect $A_{\rm CP}$ and $\gamma _{\rm CP}$. We do the analysis for three different values of the EoS, $w=0,1/3$, and $1$. For each value of $w$, we fix the end of the nonstandard epoch at $T_1=10$ MeV, which corresponds to the frequency $f(T_1)=10^{-10}$ Hz. Therefore, we are now left with only two free parameters: $H_{\rm inf}$ and $n_t$, and we do a parameter scan for them. We vary $n_t \in (-3,5)$ and $ \log_{10}\big(\frac{H_{{\rm inf}}}{\rm GeV}\big) \in (2,30)$ and compute the corresponding $A_{{\rm CP}}$ and $\gamma_{{\rm CP}}$ values. Out of the resulting points in the $A_{{\rm CP}}$-$\gamma_{{\rm CP}}$ plane, only those which lie within the 2-$\sigma$ region of the NANOGrav results have been retained.
These allowed values of $A_{{\rm CP}}$ and $\gamma_{{\rm CP}}$ are shown in Fig.~\ref{fig:Hinf_nt_allowed_rd_md} for the three different values of $w$ in there different shaded contours (gray for $w=1$, orange for $w=1/3$, and blue for $w=0$).

In Fig.~\ref{fig:Hinf_nt_allowed_rd_md}, $H_{{\rm inf}}$ values excluded from the CMB observations are shown in the magenta shaded region. It is evident from Fig.s~ \ref{fig:Hinf_nt_allowed_rd_mdw13} and~\ref{fig:Hinf_nt_allowed_rd_mdw0} that there is some part of $H_{{\rm inf}}-n_t$ parameter space which is allowed by the CMB bound for both  $w=1/3$ and $w=0$ cases. However, the parameter space required to explain the {\rm NANOGrav}~ observation for $w=1$ is fully excluded by the bound on $H_{{\rm inf}}$ from CMB (see~Fig.~\ref{fig:Hinf_nt_allowed_rd_mdw1}). Therefore, looking at the allowed values of $n_t$ for $w=1/3$ and $w=0$ in this figure, a blue tilted primordial tensor power spectrum is required to explain the NANOGrav observation in case of the standard cosmological history as well as in the presence of an early MD era.
The simplest inflationary models, like single field slow roll inflation, predict a slghtly red tilted tensor spectrum spectrum. However, a blue tilted power spectrum can arise in inflation models like G-inflation~\cite{Kobayashi:2010cm}, super-inflation models~\cite{Baldi:2005gk}, partcle produntion during inflation~\cite{Cook:2011hg,Mukohyama:2014gba}, inflation with non-Bunch Davies initial state~\cite{Ashoorioon:2014nta} and inflation with beyond slow roll~\cite{Gong:2014qga}.}

As GW contributes to the radiation energy density in the universe, it will contribute to the effective number of relativistic degrees of freedom ($N_{\rm eff}$) in the universe. An overly abundant GW can change the expansion rate too much in the early universe and therefore, can affect the abundance of light elements produced during the BBN~\cite{Allen:1997ad,Smith:2006nka,Boyle:2007zx,Kuroyanagi:2014nba,Ben-Dayan:2019gll}. BBN and Planck CMB observations constrain the contribution of the additional relativistic degrees of freedom, $\Delta N_{\rm eff}$, which requires
\begin{eqnarray}
\int df\,\frac{\Omega_{\rm GW}(f)h^2}{f} \leq 5.6 \times 10^{-6} \Delta N_{\rm eff}\, ,
\label{eq:bbncons}
\end{eqnarray}
where lower limit of the integration is equal to the frequency corresponding to the mode entering the horizon at the time of BBN and the upper limit is set to $10^7\,{\rm Hz}$, approximately corresponding to the comoving horizon at the time of the end of the inflation at temperature $T \approx 10^{15}\,{\rm GeV}$~\cite{Maggiore:1999vm,Kuroyanagi:2014nba,Liu:2015psa}. Planck CMB measurements and BBN observations severely constrain the $\Delta N_{\rm eff} \lesssim 0.4$~\cite{Aver:2015iza,Cooke:2017cwo,Aghanim:2018eyx,Hsyu:2020uqb,Aiola:2020azj}. We then compute the allowed region in the $H_{{\rm inf}}-n_t$ parameter space for which Eq.~(\ref{eq:bbncons}) is satisfied with the observed value of $\Delta N_{\rm eff}$. \pri{ The resulting constraint on $H_{{\rm inf}}-n_t$ parameter space for all the three $w$ under consideration are shown by the yellow dotted line in Fig.~\ref{fig:Hinf_nt_allowed_rd_md}, where the region above that line (shaded in yellow colour) is ruled out.
In addition to the bound on GWs from their contribution to the $N_{\rm eff}$, non-observation of the stochastic GW background in LIGO O1 and O2 runs can also put the constraints on the $H_{{\rm inf}}-n_t$ parameter space~\cite{LIGOScientific:2019vic}. The constraints for all three $w$ are shown by the green dashed lines in Fig.~\ref{fig:Hinf_nt_allowed_rd_md}, where the region above the lines (green shaded region) are excluded.  To calculate the bound from LIGO, the primordial power spectrum is assumed to retain the power law form until the LIGO frequencies. It is evident from Fig.~\ref{fig:Hinf_nt_allowed_rd_mdw13} that the entire region of $H_{{\rm inf}}-n_t$ parameter space allowed by the NANOGrav and CMB bounds in the standard RD history is ruled out when the bound on $\Delta N_{\rm eff}$ from BBN and Planck CMB observations and the LIGO bound are taken into account.  However, there is a finite region left in the $H_{{\rm inf}}-n_t$  parameter space which can explain the NANOGrav result and simultaneously evades the constraints from LIGO, BBN, and CMB in the case of nonstandard thermal history with an early MD epoch (see Fig.~ \ref{fig:Hinf_nt_allowed_rd_mdw0}).  The $\gamma_{{\rm CP}}-A_{{\rm CP}}$ parameter space corresponding to the $H_{{\rm inf}}-n_t$ parameter space allowed by the {\rm NANOGrav}~, CMB, BBN, and LIGO observation in a nonstandard cosmological history with a pre-BBN MD era is shown by the blue shaded region in Fig.~\ref{fig:Hinf_nt_allowed_rd_mdw0acp}). }

\section{Induced primordial gravitational waves from curvature perturbations}
\label{GW2sec}
Scalar and tensor perturbations do not evolve independently at the second order of perturbation theory. \suk{In certain single field and multi-field models of inflation, the scalar perturbations can become large at small length scales owing to various mechanisms in the inflationary dynamics~\cite{Garcia-Bellido:2017mdw,Ballesteros:2017fsr,Ballesteros:2020qam,Bhaumik:2019tvl,Germani:2017bcs,Mishra:2019pzq,Clesse:2015wea,Braglia:2020eai,Braglia:2020taf,Fumagalli:2020adf,Palma:2020ejf,Ozsoy:2020kat,Aldabergenov:2020yok}. The resulting GWs induced by the large scalar perturbations can be large enough to be relevant in view of the (proposed) sensitivities of the current and future GW detectors. Such large scalar fluctuations can also lead to abundant formation of PBHs. Therefore, it is interesting to check the implications for PBH abundance assuming the observed {\rm NANOGrav}~ signal to be second order GW induced by the first order scalar perturbations. In this section, we analyse such implications for the standard and nonstandard post-inflationary evolutions for three different profiles of the small scale growth of the scalar power spectra. The difference between the analysis in this section and the previous one is that here, the primordial scalar perturbations have large amplitude at small scales, whereas, in the previous section the small scale primordial tensor perturbations (first order) were large. This difference directly relates to different features in the inflationary dynamics that lead to growth in either the scalar or the tensor power spectra.}

Evolution of the second order tensor perturbation is given as
\begin{equation}
\label{equ:hhh}
h_k'' + 2 {\cal H} h_k' + k^2 h_k=  {\cal S_{\mathbf{k}}}(\eta)\, ,
\end{equation}
where ${\cal S}(\mathbf{k}, \eta)$ is the source term which depends on the first order scalar perturbations. The scalar induced second order tensor power spectrum is~\cite{Kohri:2018awv},
\begin{align}\label{eq:pgamma}
\overline{\mathcal{P}_h(k,\eta)}=2\int_0^\infty dv \int_{|1-v|}^{1+v} du\left[\frac{4v^2-\left(1+v^2-u^2\right)^2}{4uv}\right]^2 P_{ \zeta}(kv)P_{\zeta}(ku)\overline{I^2}(v,u,x)\,,
\end{align}
where the bar denotes the oscillation average.
Here $v\equiv q/k$, $u\equiv|\mathbf{k}-\mathbf{q}|/k$, $P_{ \zeta}$ is the primordial curvature power spectrum, and $x=k\eta$. The integration kernel $I(u,v,x)$ contains the source information (complete expression can be found in ref~\cite{Espinosa:2018eve,Kohri:2018awv}) which can be written with a simple analytical expression for very late time, i.e. $x \gg 1$. For the modes entering the horizon during RD era,
\begin{align}
\overline{I_{RD}^2}(x\gg 1,u,v)=\frac{9(u^2+v^2-3)^2}{32u^6v^6x^2}&\Bigg\{ {\pi}^2(u^2+v^2-3)^2\Theta(u+v-\sqrt{3}) \nonumber\\&+  \bigg(-4uv + (u^2+v^2-3) \ln {\frac{3-(u+v)^2}{3-(u-v)^2} } \bigg)^2\Bigg\}\,.
\label{IRD}
\end{align}
For a general EoS ($0 < w \leq 1$), the expression for $I(u,v,x)$ is given as~\cite{Domenech:2019quo}
\begin{align}\label{eq:kernel3}
I(x\gg1,u,v)&=2^{\beta}\frac{3\sqrt{2}}{w\pi\alpha^3}\frac{1+w}{1+3w}\Gamma^2[\beta+2]{\left(uvx\right)^{-\beta-1/2}}\nonumber\\&\times\left\{\frac{\pi}{2}\sin\left(x-\frac{\beta\pi}{2}-\frac{\pi}{4}\right){I}_J(u,v,w)+\cos\left(x-\frac{\beta\pi}{2}-\frac{\pi}{4}\right){I}_Y(u,v,w)\right\}\, ,
\end{align}
where $I_J(u,v,w)$ and $I_Y(u,v,w)$ are expressed in terms of the associated Legendre polynomials and the Legendre polynomials on the cut (complete expression can be found in ref~\cite{Domenech:2019quo}). In general, it is numerically challenging to calculate the second order tensor power spectrum for an arbitrary value of EoS ($0 < w \leq 1$) due to the complicated dependence of $I_J$ and $I_Y$ on $w$. Hence, we chose two different values to interpret the dynamics in nonRD epochs: (i) $w=1$, where the energy density falls faster than radiation and (ii) $w=\frac{1}{9}$, where it falls slower. Rest of the analysis is carried out for these two nonstandard epochs in addition to the standard RD epoch. EoS $w=1$ can arise in models of quintessential inflation~\cite{Peebles:1998qn}, where, after the end of inflation, the inflaton fast rolls (kinetic energy dominated (KD) universe) to reach the potential for quintessence dark energy. An epoch with $w=1/9$ is a nearly dust-like EoS, which can arise during very slow reheating process of the inflaton or the moduli such that the effective EoS can lie between MD and RD for some considerable time. The particular choice of $w=1/9$ is because the expressions for $I_J$ and $I_Y$ are much simplified in this case and can be solved analytically for $x\gg 1$. The nonstandard epochs terminate at $T_1\geq T_{\rm BBN}$ and we choose $T_1=10$ MeV.

The integral $I_{{\rm KD}}$ for $w=1$ is estimated as:
\begin{align}\label{IKD}
\overline{I_{KD}^2}(x\gg1,u,v)&=\frac{9}{16\pi u^4v^4x}\left\{\frac{\left(3(u^2+v^2-1)^2-4u^2v^2\right)^2}{{4u^2v^2-(u^2+v^2-1)^2}}+9\left(u^2+v^2-1\right)^2\right\}\,.
\end{align}
Similarly, the expression for $I_{1/9}$ can be derived for $w=\frac{1}{9}$ and given as
\begin{align}
\label{I19D}
\overline{I^2_{1/9}}(x\gg1,u,v)&=\frac{25}{16\pi u^8v^8{x}^{3}}\Bigg\{\left(u^2+v^2-9\right)^2\left[5\left(u^2+v^2-9\right)^2-6u^2v^2\right]^2\nonumber\\&
+{\left|4u^2v^2-(u^2+v^2-9)^2\right|}\left(4u^2v^2+5(u^2+v^2-9)^2\right)^2\nonumber\\&+
2\left(u^2+v^2-9\right)\left[5\left(u^2+v^2-9\right)^2-6u^2v^2\right]\sqrt{\left|4u^2v^2-(u^2+v^2-9)^2\right|}\nonumber\\&
\times\left(4u^2v^2+5(u^2+v^2-9)^2\right)\Theta(\sqrt{9}-u-v)\Bigg\}\,.
\end{align}
Now, using the expressions for the integration kernels in Eq.s~\eqref{IRD},~\eqref{IKD} and~\eqref{I19D} with the Eq.s~\eqref{originalGWPh} and~\eqref{eq:pgamma}, we can determine the {$\Omega_{\rm GW}^{(2)}$}~ for standard RD and the two nonstandard EoS. For nonstandard EoS, the respective $I_{w}$ contribute to {$\Omega_{\rm GW}^{(2)}$}~ until the temperature $T_1$ (corresponding to the comoving time $\eta_1$) is reached, after which, standard evolution with RD follows\footnote{Here, we do not explicitly solve for the continuity of {$\Omega_{\rm GW}^{(2)}$}~ at $\eta _1$, since the scales corresponding to the transition are outside the {\rm NANOGrav}~ relevant scales. Moreover, we consider the transition from the nonstandard to RD epoch to be instantaneous. In a more generic analysis with a gradual transition between these two epochs, it is important to account for the continuous transition of {$\Omega_{\rm GW}^{(2)}$}~ at $\eta _1$.}. Thus, for a general $w$ dominated epoch~\cite{Domenech:2019quo},
\begin{equation}
 \Omega_{{\rm GW}}^{(2)}(k,\eta _0)= 0.387 ~\Omega _{{\rm rad},0} \bigg(\frac{g_s^4 g_*^{-3}}{106.75}\bigg)^{-1/3} \Omega_{{\rm GW}}^{(2)} (k,\eta _1),
\label{omg2fullnow}
\end{equation}
where
\begin{equation}
 \Omega_{{\rm GW}}^{(2)} (k,\eta _1)=\frac{(1+3w)^2}{48}(k\eta _1)^2 \overline{\mathcal{P}_h(k,\eta _1)}. \label{omg2trans}
\end{equation}
Therefore, a general $I_{w}$ contributes to Eq.~\eqref{omg2trans} and the full second order GW can be calculated then using Eq.~\eqref{omg2fullnow} once the primordial power spectrum $P_{\zeta}(k)$ is specified.

\subsection{PBH formation}
\label{secPBH}
If the amplitude of the primordial fluctuatiions is such that when the modes re-enter the horizon at the post-inflationary epoch, the density fluctuations $\delta$ are larger than the critical density for collapse ($\delta_c$), then PBH can be produced with mass $M=\epsilon M_H$, where $M_H$ is the horizon mass at collapse. $\epsilon$ signifies the efficiency of collapse and is typically of order $1$; here we consider $\epsilon = 0.33$. The PBH formation process depends on the background EoS through the dependence of $\delta_c$ on $w$ and the evolution of PBH abundance in this epoch depends on the modified expansion rate. Taking the numerically evaluated expression in~\cite{Harada:2013epa}:
\begin{equation}
\delta _c(w)=\frac{3(1+w)}{(5+3w)}\sin ^2 \bigg(\frac{\pi \sqrt{w}}{(1+3w)}\bigg). \label{deltac_w}
\end{equation}
Considering the Press-Schechter formalism for gravitational collapse, present abundance of PBH in an interval of mass $M$ to $M+dM$ for general $w$ is given as~\cite{Bhattacharya:2019bvk}
\begin{equation}
\psi_w(M)=\frac{\epsilon}{T_{\rm eq}}\bigg(\frac{g_s(T_1)}{g_s(T_{\rm eq})}\bigg)^{\frac{1}{3}}\bigg(\frac{\Omega_m h^2}{\Omega_c h^2}\bigg)\bigg(\frac{90M_P^2}{\pi ^2g_*(T_1)}\bigg)^{\frac{w}{1+w}}(4\pi \epsilon M_P^2)^{\frac{2w}{1+w}}T_1^{\frac{1-3w}{1+w}}\frac{\beta (M)}{M^{\frac{3w+1}{1+w}}},\label{psiM}
\end{equation}
where $\epsilon$ is the fraction of the total horizon mass that contributes to the PBH mass~\cite{Carr:1975qj} and $M_P$ is the reduced Planck mass. The PBH mass fraction is given as:
\begin{equation}
\beta(M)={\rm erfc} \bigg[\frac{\delta_c(w)}{\sqrt{2\sigma_{\delta}^2}}\bigg], \label{beta_w}
\end{equation}
where $\sigma_{\delta}^2$ is the variance of the density power spectrum and calculated as:
\begin{equation}
\sigma_{\delta}^2=\frac{4(1+w)^2}{(5+3w)^2}\int \frac{dk}{k}(kR)^4 W^2(k,R)P_{\zeta}(k).
\label{vardel}
\end{equation}
We choose the window function $W^2(k,R)=\exp(\frac{-k^2R^2}{4})$ to smooth the perturbations on the comoving scale $R$ at formation. The mass $M$ of the PBH produced is related to the comoving wavenumber $k$ at formation for general EoS as:
\begin{equation}
M(k)=4\pi\epsilon {M_P}^2\bigg(\frac{\pi ^2 g_*^{\rm eq}}{45M_P^2}\bigg)^{\frac{1}{3w+1}}\bigg(\frac{g_s^{\rm eq}}{g_s(T_1)}\bigg)^{\frac{3w-1}{3(3w+1)}}(a_{\rm eq}T_{\rm eq})^{\frac{3(1+w)}{3w+1}}T_1^{-\frac{3w-1}{3w+1}}k^{-\frac{3(1+w)}{3w+1}}. \label{Mk_w}
\end{equation}
$\psi_w(M)$ can be calculated from Eq.~\eqref{psiM} using Eq.s~\eqref{deltac_w},~\eqref{beta_w},~\eqref{vardel} and~\eqref{Mk_w}. The ratio of the PBH abundance for a general $w$ dominated epoch and that for a RD epoch is:
\begin{equation}
g_{\psi}(M,w)\equiv \frac{\psi_w(M)}{\psi_{\rm RD}(M)}=\frac{\beta (M)}{\beta ^{\rm RD}(M)}\frac{g_s(T)^{w-1/3}}{g_s(T_1)^{w-1}g_{s}(T_{\rm eq})^{2/3}}\bigg(\frac{T}{T_1}\bigg)^{3w-1}.\label{gain01}
\end{equation}
The dependence of $\delta_c$ on $w$ is such that $\delta_c^{(w=1)},\delta_c^{(w=1/9)} < \delta_c^{(w=1/3)}$ and therefore, the gain $g_{\psi}(M,w)>1$ for both the nonstandard epochs. This means that to form same total abundance $f_{\rm PBH}=\int dM \psi_w(M)$, the two nonstandard epochs will require a lower value of $A_s$ than that required in RD. In this work, we check if the low value of $A_s$ required for abundant PBH formation in each nonstandard epoch can satisfy the {\rm NANOGrav}~ bounds.
\subsection{Primordial scalar power spectrum}
\suk{An enhanced $P_{\zeta}(k)$ at scales smaller than that observed in CMB surveys, can be motivated in a variety of inflation models. In single-field inflation models, presence of an inflection point or a bump can lead to ultra slow-roll of the inflaton, which results in enhancement of $P_{\zeta}(k)$~\cite{Garcia-Bellido:2017mdw,Ballesteros:2017fsr,Ballesteros:2020qam,Bhaumik:2019tvl,Germani:2017bcs,Mishra:2019pzq}. For multi-field models of inflation, the inflaton coupled to a secondary field (e.g. mild waterfall phase in hybrid inflation~\cite{Clesse:2015wea} or noncanonical coupling~\cite{Braglia:2020eai,Braglia:2020taf}) or rapid turns in field space~\cite{Fumagalli:2020adf,Palma:2020ejf} can induce large $P_{\zeta}(k)$.
However, in the rest of this paper, we proceed with an agnostic approach towards the exact inflationary model that gives rise to an enhanced $P_{\zeta}(k)$; rather, we consider three well-motivated functional forms for $P_{\zeta}(k)$ and analyse the PBH abundance and relevant second order GW spectrum for each case. Throughout the analysis of this paper, we consider the scalar fluctuations to be Gaussian, and any nontrivial non-Gaussianity can affect the results for PBH abundance and SGWB spectra and therefore alter the implications derived here.}
\subsubsection{Gaussian power spectrum}
\label{GPP}
A Gaussian primordial power spectrum is motivated in many single field and multi-field inflation models discussed above and is given as:
\begin{equation}
P_{\zeta}(k)=A_s\times \exp\bigg[-\frac{(\log({k/k_*}))^2}{2\sigma_p^2}\bigg]\,.
\label{PKG}
\end{equation}
The GW spectrum, normalised with $A_s^2$ is given in Fig.~\ref{OmgA2}.
\begin{figure}[!htbp]
\centering{
\includegraphics[width=0.6\textwidth]{omgwbyA2_final.pdf}
\caption{\small{The normalised {$\Omega_{\rm GW}^{(2)}$}~ for the three $w$ cases are shown as a function of $k/k_*$.}}
\label{OmgA2}}
\end{figure}
Evidently, for $w=1$, there is an enhancement in the normalised GW spectra and for $w=1/9$, it is diminished. Moreover, the peak of the spectrum shifts towards $k>k_*$ for $w=1$ and towards $k<k_*$ for $w=1/9$. Now, the {\rm NANOGrav}~ analysis fits the observed stochastic spectrum in the frequency range $f_0\in (2.5\times10^{-9},1.2\times10^{-8})$ Hz. Thus, if we assume that the {$\Omega_{\rm GW}^{(2)}$}~ peak for each $w$ occurs at the scale (wavenumber $k_0=2\pi f_0$) observed by the {\rm NANOGrav}~, then for $w=1/3$, $k_0/k_*\simeq 1$; for $w=1$, $k_0/k_*>1$, and for $w=1/9$, $k_0/k_*<1$. However, in that case, the amplitudes required for each case to fit the {\rm NANOGrav}~ result will be different: $A_s(w=1/3)=6\times 10^{-3}$, $A_s(w=1)=5.5\times 10^{-4}$, and $A_s(w=1/9)=5.1\times 10^{-3}$ to have {$\Omega_{\rm GW}^{(2)}$}~$_{\rm ,peak}\simeq 10^{-9}$.
\begin{figure}[!htbp]
\centering{
\includegraphics[width=0.6\textwidth]{PBH_psi.pdf}
\caption{\small{PBH mass functions for different $w$ epochs are shown for a typical $A_s=0.007$. Colour specifications are same as Fig.~\ref{OmgA2}.}}
\label{psiPBH}}
\end{figure}

The PBH mass function $\psi_w(M)$ using this $P_{\zeta}(k)$ for $w=1/3$, $1$, and $1/9$ are evaluated and given in Fig.~\ref{psiPBH} for $A_s=0.007$. This figure shows an enhancement in PBH abundance for the nonRD epochs and peak shift with respect to peak position in the RD. These peak shifts are correlated with the shifts in GW in Fig.~\ref{OmgA2}.
\begin{figure}[!htbp]
\centering{
\includegraphics[width=0.75\textwidth]{dotsGandC_rightprefacnT1.pdf}
\caption{\small{$A_{\rm CP}$ and $\gamma _{\rm CP}$ fitted for $f_{\rm PBH}=0.1$ for each case of a Gaussian primordial spectrum: $\sigma_p=1$, $2$, and $3$ are plotted with circle, square and triangles respectively;  for $w=1/3$ in red,  $w=1$ in blue and  $w=1/9$ in green. Points for $f_{\rm PBH}=0.1$ for a constant primordial power spectrum for these three $w$ are plotted as pink stars. This constant power spectrum generates PBH masses in the range: $0.1<M/M_{\odot}<5$. The dashed gray line signifies the open edge of the $2\sigma $ contour in the {\rm NANOGrav}~ analysis (orange contour in Fig.1 of~\cite{Arzoumanian:2020vkk}).}}
\label{dotplot}}
\end{figure}
For the curvature power spectrum given in Eq.~\eqref{PKG} with $\sigma_p=$1, 2, and 3, we evaluate {$\Omega_{\rm GW}^{(2)}$}~ with $k_*=3.6\times 10^6$ Mpc$^{-1}$ which corresponds to $f_*=5.5$ nHz, and fit the results in the frequency range of the {\rm NANOGrav}~ for each of the three $w$ epochs. For each case, we chose $A_s$ such that $f_{\rm PBH}=0.1$. The resulting $A_{\rm CP}$ and $\gamma _{CP}$ are plotted with the {\rm NANOGrav}~ $1\sigma$ and $2\sigma$ contours in Fig.~\ref{dotplot}.  For RD (red circle, square, and triangle), the resulting three points  for different $\sigma_p$ are close together inside the contour, therefore, abundant PBH production is consistent with the {\rm NANOGrav}~ observation. However, for the other two cases ($w=1$ in blue and $w=1/9$ in green), points for different $\sigma_p$ are far apart with the value closest to the contours given by $\sigma_p=2$ for both the cases. Therefore, in the rest of our analysis, we use $\sigma_p=2$ to check the consistency of abundant PBH formation in these nonRD epochs with the {\rm NANOGrav}~ result. In this figure, we also show the points corresponding to a constant power spectrum of the following form:
\begin{equation}
P_{\zeta}(k)=A_s\times \Theta (k-k_l)\times \Theta (k_s-k),
\label{constpow}
\end{equation}
where $k_l$ and $k_s$ correspond to the scales between which the power spectrum is constant. We have chosen $k_l$ and $k_s$ values for each $w$ case such that, the resulting PBHs from~\eqref{constpow} are always in the mass range $M/M_{\odot}\in (0.1,5)$. We find that for the constant power spectrum, PBH production with $f_{\rm PBH}=0.1$ (stars in Fig.~\ref{dotplot}) is not consistent with the {\rm NANOGrav}~ result in these nonstandard epochs under consideration.
\begin{table}[!htbp]
\caption{Range of amplitude $A_s$ chosen for the three cases with different EoS used in Fig.~\ref{regcont}.}
\begin{center}
\label{TT1}
\resizebox{0.4\textwidth}{!}{
\begin{tabular}{|c|c|c|c|}
\hline
 & $w=1/3$ & $w=1$ & $w=1/9$ \\
\hline
\hline
$A_s^{\rm max} $ & 0.015 & 0.007 & 0.0082\\
\hline
$A_s^{\rm min}$ & 0.002 & 0.001 & 0.003\\
\hline
\end{tabular}}
\end{center}
\end{table}

Fixing $\sigma_p=2$, we vary $A_s$ (ranges given in Table~\ref{TT1}) and vary $k_*$ in the range ($2\times 10^6$ Mpc$^{-1}$ - $7 \times 10^6$ Mpc$^{-1}$) for each $w$ scenario. The corresponding variations in $A_{\rm CP}$-$\gamma_{\rm CP}$ space are given by magenta ($w=1/3$), cyan ($w=1$), and green ($w=1/9$) shaded regions in Fig.~\ref{regcont}. We show the contours for PBH abundance in dashed lines inside these contours with $f_{\rm PBH}=0.1$, $10^{-5}$, and $10^{-10}$ from top to bottom. It is evident from Fig.~\ref{regcont} that $10\%$ PBH as dark matter is possible for RD. For $w=1$, All these dashed contours remain outside of the {\rm NANOGrav}~ $2\sigma$ contour; the most hopeful case in this $k_*$ range is represented by the cyan dotted contour at the bottom, which represents $f_{\rm PBH}\simeq 10^{-22}$. Therefore, abundant PBH production in $w=1$ epoch is not consistent with the {\rm NANOGrav}~ observation. Increasing $k_*$ will lead to the {$\Omega_{\rm GW}^{(2)}$}~ peak being closer to the {\rm NANOGrav}~ frequency range, and the $f_{\rm PBH}$ contours come down to lower $A_{\rm CP}$ values. However, the slope of {$\Omega_{\rm GW}^{(2)}$}~ is such that $f_{\rm PBH}=0.1$ contour comes within the {\rm NANOGrav}~ $A_{\rm CP}$ range only when $\gamma_{\rm CP}$ is very small. Moving towards the lower values of $k_*$ does not help since power at $k_0$ goes further down and will also influence the power spectrum at CMB scales for a broad primordial spectrum.
\begin{figure}[!htbp]
\centering{
\includegraphics[width=0.75\textwidth]{Region_contour01_final.pdf}
\caption{\small{$A_{\rm CP}$ and $\gamma _{\rm CP}$ fitted for $A_s$ varied in different ranges for $w=1/3$ (magenta shaded), $w=1$ (cyan shaded), and $w=1/9$ (green shaded) with variation in the pivot scale in the range $k_*=2\times 10^6$ Mpc$^{-1}$ to $k_*=7\times 10^6$ Mpc$^{-1}$ for the Gaussian primordial spectrum with $\sigma_p=2$. The dashed lines inside each shaded region are contours for $f_{\rm PBH}=0.1$, $10^{-5}$ and $10^{-10}$ from top to bottom. The additional dotted cyan line in the $w=1$ case just enters the 2$\sigma$ contour of the {\rm NANOGrav}~ observation, but does not motivate abundant PBH (here $f_{\rm PBH}\simeq 10^{-22}$). For the $w=1/9$ case, the upper edge of the green shaded region motivates $f_{\rm PBH}=1$, but the observations constrain $f_{\rm PBH}$ to a much lower value.}}
\label{regcont}}
\end{figure}
For $w=1/9$, the $f_{\rm PBH}=0.1$ contour enters the {\rm NANOGrav}~ contour only for $k_*\simeq 7 \times 10^6$ Mpc$^{-1}$, and therefore it is possible to achieve abundant PBH production for the $w=1/9$ case for $k_*$ larger than the {\rm NANOGrav}~ probed range.
Therefore, we proceed with the case $w=1/9$ now and vary $k_*$ in the range ($2\times 10^7$ Mpc$^{-1}$ - $7 \times 10^7$ Mpc$^{-1}$) for a variation in $A_s$ in the range $0.005-0.0082$ for $w=1/9$. For the RD case, $A_s$ is varied in the same $k_*$ range as before, but for a smaller range in $A_s$ between $0.008-0.015$. The resulting variation in $A_{\rm CP}$-$\gamma_{\rm CP}$ space are shown in Fig.~\ref{regcont2} in green for $w=1/9$ and in magenta for $w=1/3$. The relevant PBH contours for $f_{\rm PBH}=0.1$ are shown in dashed lines. We conclude that the PBH and the second order GW production in a nonstandard epoch with EoS $w=1/9$ is consistent with the {\rm NANOGrav}~ result when $k_*$ is away from the {\rm NANOGrav}~ probed range. The gain in the $w=1/9$ over RD epoch is that for this nonstandard epoch, abundant PBH and observed GW can both be explained with a lower value of $A_s$: $A_s(w=1/9)\vert _{f_{\rm PBH}=0.1}\simeq 0.007$ and $A_s(w=1/3)\vert _{f_{\rm PBH}=0.1}\simeq 0.01$.
\begin{figure}[!htbp]
\centering{
\includegraphics[width=0.75\textwidth]{Region_contour02_final.pdf}
\caption{\small{For a Gaussian primordial power spectrum, $A_{\rm CP}$ and $\gamma _{\rm CP}$ fitted for $A_s$ varied in different ranges for $w=1/3$ (magenta shaded) and $w=1/9$ (green shaded) with variation in the pivot scale in the range $k_*=2\times 10^6$ Mpc$^{-1}$ to $k_*=7\times 10^6$ Mpc$^{-1}$ for RD and $k_*=2\times 10^7$ Mpc$^{-1}$ to $k_*=7\times 10^7$ Mpc$^{-1}$ for $w=1/9$. The dashed line inside each shaded region is the $f_{\rm PBH}=0.1$ contour. The value of $\sigma_p$ is fixed to be 2.}}
\label{regcont2}}
\end{figure}
The Gaussian power spectrum with $\sigma_p=2$ has a broad peak, therefore we quote a range of the PBH mass near the peaks for each case for the lowest and highest $k_*$ under consideration in Table~\ref{TT2}.
\begin{table}[!htbp]
\caption{Ranges of PBH mass $M$ for each EoS}
\begin{center}
\label{TT2}
\resizebox{0.85\textwidth}{!}{
\begin{tabular}{|c|c|c|c|}
\hline
 & $w=1/3$ & $w=1$ & $w=1/9$ \\
\hline
\hline
Range of $k_*$ in Mpc$^{-1}$ & $2\times 10^6$ - $7 \times 10^6$ & $2\times 10^6$ - $7 \times 10^6$  & $2\times 10^7$ - $7 \times 10^7$\\
\hline
Range of $M/M_{\odot}$ at $k_{*,{\rm min}}$ & 0.2-2 & 1-10 & $10^{-4}$ - $5\times 10^{-3}$\\
\hline
Range of $M/M_{\odot}$ at $k_{*,{\rm max}}$ & 0.01 - 0.33 & 0.08-2 & $3\times 10^{-6}$ - $3\times 10^{-4}$\\
\hline
\end{tabular}}
\end{center}
\end{table}
\subsubsection{Lognormal power spectrum with an exponential cut off}
\label{PLN}
\suk{The lognormal power spectrum with an exponential cutoff is motivated in~\cite{Braglia:2020eai,Vaskonen:2020lbd} and is given as
\begin{equation}
P_{\zeta}(k)=A_s\times \exp\bigg[\beta \bigg(1-\frac{k}{k_*}+\log \bigg(\frac{k}{k_*}\bigg) \bigg)-\alpha \log ^2 \big(\frac{k}{k_*}\big) \bigg]\,,
\label{PKLN}
\end{equation}
where $\alpha=0.17$ and $\beta=0.62$. The second order GW spectra normalized by $A_s^2$ for this power spectrum are shown in Fig.~\ref{OmgA2_LN} for the standard RD epoch and other two nonstandard post-inflationary epochs with $w=1$ and $w=1/9$.
\begin{figure}[!htbp]
\centering{
\includegraphics[width=0.6\textwidth]{Omg_A2_LN.pdf}
\caption{\small{The normalised {$\Omega_{\rm GW}^{(2)}$}~ for the three values of $w$ are shown as a function of $k/k_*$ using the lognormal power spectrum with an exponential cutoff.}}
\label{OmgA2_LN}}
\end{figure}
The profiles of the normalised \OMG2in Fig.~\ref{OmgA2_LN} show that the peaks occur at $k=k_*$ for RD, $k<k_*$ for $w=1/9$, and $k\simeq 1.38k_*$ for $w=1$. Assuming that the location of the peaks in the {$\Omega_{\rm GW}^{(2)}$}~ spectra are within the {\rm NANOGrav}~ frequency range, the amplitude of $P_{\zeta}(k)$ required for each $w$ to explain the observed GW amplitude of $10^{-9}$ are $A_s (w=1/3)=6.6\times 10^{-3}$, $A_s (w=1)=9.5\times 10^{-4}$, and $A_s (w=1/9)=0.0146$.

The PBH mass function for this $P_{\zeta}(k)$ is plotted in Fig.~\ref{psiPBH_LN} for the three $w$ values, which show that each of them corresponds to a broad range of PBH mass. Here, with the given choices of $k_*=k_0$ and $A_s$, the PBH masses relevant to the binary merger events observed in LIGO/Virgo are produced with the most significant abundance only for $w=1$, whereas the $w=1/9$ case permits a considerable abundance of the planetary mass PBHs.
\begin{figure}[!htbp]
\centering{
\includegraphics[width=0.6\textwidth]{mpsim_LN.pdf}
\caption{\small{PBH mass functions for different $w$ epochs are shown for a typical primordial amplitude $A_s=0.007$ using the power spectrum given in Eq.~\ref{PKLN}. Colour specifications are same as Fig.~\ref{OmgA2_LN}.}}
\label{psiPBH_LN}}
\end{figure}

A similar analysis as in Fig.~\ref{regcont} in Sec.~\ref{GPP} is carried out here for the power spectrum in Eq.~\ref{PKLN} where $k_*$ is varied in the range $(10^6 - 10^7)$ Mpc$^{-1}$ and $A_s$ is varied in the range $(0.002 - 0.02)$. The resulting scanned parameter space in the $A_{\rm CP} - \gamma_{{\rm CP}}$ plane is shown in Fig.~\ref{LNcont} for $w=1/3$ (magenta), $w=1$ (cyan), and $w=1/9$ (green), where the dark shaded region in each case corresponds to the exact frequency range where the {\rm NANOGrav}~ analysis fits the observed spectrum as a common power law. The coloured dashed contours in each case signify PBH abundance of $0.1$ (upper) and $10^{-5}$ (lower). Embedding the $1\sigma$ and $2\sigma$ contours for $A_{\rm CP} - \gamma_{{\rm CP}}$ from the {\rm NANOGrav}~ observation in this plot (black contours), it is evident that for RD, $10\%$ PBH abundance for $P_{\zeta}(k)$ given in Eq.~\ref{PKLN} is consistent with the {\rm NANOGrav}~ 12.5 year data at $2\sigma$ confidence level. The nonstandard scenario with $w=1/9$ is also consistent with the {\rm NANOGrav}~ data where PBH is produced with $\gtrsim 10\%$ abundance, but only for $k_*\gtrsim 5\times 10^6$ Mpc$^{-1}$. On the other hand, the second order GW produced in the nonstandard $w=1$ epoch is consistent with the {\rm NANOGrav}~ observation only when the PBH abundance is negligible.
\begin{figure}[!htbp]
\centering{
\includegraphics[width=0.75\textwidth]{LNshade_withpbhcont.pdf}
\caption{\small{For a lognormal primordial power spectrum with an exponential cut off, $A_{\rm CP}$ and $\gamma _{\rm CP}$ fitted for $A_s$ varied in the range $0.002 - 0.02$ for $w=1/3$ (magenta shaded), $w=1$ (cyan shaded), and $w=1/9$ (green shaded). For dark shaded regions, the pivot wavenumber runs from $k_*=1.615\times 10^6$ Mpc$^{-1}$ to $k_*=7.756\times 10^6$ Mpc$^{-1}$ from right to left, where {\rm NANOGrav}~ fits their observation to a common power law power spectrum. The light shaded regions are shown for the pivot wavenumber $k_*=10^6$ Mpc$^{-1}$ to $k_*=10^7$ Mpc$^{-1}$ from right to left. The dashed contours are for the PBH abundance of $0.1$ (above) and $10^{-5}$ (below) for each case. Black solid contours represent the $1\sigma$ and $2\sigma$ confidence levels for the parameters $A_{\rm CP}$ and $\gamma _{\rm CP}$ observed in the {\rm NANOGrav}~ .}}
\label{LNcont}}
\end{figure}}

\subsubsection{Broken power law power spectrum}
\label{PBPL}
\suk{The primordial power spectrum can also have the form of a broken power law~\cite{Ballesteros:2017fsr,Byrnes:2018txb,Vaskonen:2020lbd} given by
\begin{equation}
P_{\zeta}(k)=A_s\times \frac{\alpha +\beta}{\beta \bigg(\frac{k}{k_*}\bigg)^{-\alpha}+\alpha \bigg(\frac{k}{k_*}\bigg)^{\beta}}\, ,
\label{PKBPL}
\end{equation}
where $\alpha=4$ and $\beta=0.6$. The second order GW spectra normalized by $A_s^2$ for this power spectrum are shown in Fig.~\ref{OmgA2_BPL} for the epochs $w=1/3$, $1$, and $1/9$. Unlike the power spectra discussed in Sec.~\ref{GPP} and Sec.~\ref{PLN}, this form of $P_{\zeta}(k)$ has a sharp peak, which is reflected in the wedge-like shapes of the normalized {$\Omega_{\rm GW}^{(2)}$}~ in Fig.~\ref{OmgA2_BPL}. Moreover, the difference in slopes of {$\Omega_{\rm GW}^{(2)}$}~ among the three $w$ cases clearly manifest the power-law dependence of the integration kernel $\overline{I^2}(v,u,x)$ in Eq.~\ref{eq:pgamma} on $w$.  The peak positions in this figure are $k=k_*$ for $w=1/3$, $k\simeq 3k_*$ for $w=1$, and $k\simeq 0.5k_*$ for $w=1/9$ respectively. Assuming that for each $w$, the peaks occur in the {\rm NANOGrav}~ observed frequency range, the peak amplitudes required to explain the {\rm NANOGrav}~ observation with {$\Omega_{\rm GW}^{(2)}$}~ $\sim 10^{-9}$ are $A_s=6.9\times 10^{-3}$, $8.7\times 10^{-4}$, and $0.0187$ for $w=1/3$, $1$, and $1/9$ respectively.

\begin{figure}[!htbp]
\centering{
\includegraphics[width=0.6\textwidth]{Omg_A2_BPL.pdf}
\caption{\small{The normalised {$\Omega_{\rm GW}^{(2)}$}~ for the three values of $w$ are shown as a function of $k/k_*$ using the broken power law power spectrum.}}
\label{OmgA2_BPL}}
\end{figure}
\begin{figure}[!htbp]
	\centering{
		\includegraphics[width=0.6\textwidth]{mpsim_BPL.pdf}
		\caption{\small{PBH mass functions for different $w$ epochs are shown for a typical primordial amplitude $A_s=0.007$ using the power spectrum given in Eq.~\ref{PKBPL}. Colour specifications are same as Fig.~\ref{OmgA2_BPL}.}}
		\label{psiPBH_BPL}}
\end{figure}
The corresponding PBH mass function is shown in Fig.~\ref{psiPBH_BPL} for all the epochs under consideration. It is evident from Fig.~\ref{psiPBH_BPL} that PBHs are produced in different mass ranges for each case. Similar to the results in Sec.~\ref{GPP} and Sec.~\ref{PLN}, PBHs in the LIGO/Virgo observed mass range are produced with the large abundance only for $w=1$, whereas the planetary mass PBHs are produced in a significant number for $w=1/9$.

The implications of the broken power law power spectrum in Eq.~\ref{PKBPL} for {$\Omega_{\rm GW}^{(2)}$}~ and PBH abundance are shown in the $A_{\rm CP} - \gamma_{{\rm CP}}$ plane in Fig.~\ref{BPLcont} with the color and line specifications similar to Fig.s~\ref{regcont2} and~\ref{LNcont}. Varying $k_*$ in the range $(10^6 - 10^7)$ Mpc$^{-1}$ and $A_s$  in the range $(0.002 - 0.02)$ show that if the {\rm NANOGrav}~ observation is assumed to be due to the second order GW production in the nonstandard $w=1$ epoch, then PBH is produced negligibly. For the standard RD epoch, PBH with $10\%$ abundance can be produced here while the corresponding {$\Omega_{\rm GW}^{(2)}$}~ is consistent with {\rm NANOGrav}~ result within $1\sigma$ confidence limit. However, for the nonstandard dustlike case with $w=1/9$, even $100\%$ PBH abundance (although observationally not allowed) can be reached with an amplitude $A_s$ lower than that required to match predicted {$\Omega_{\rm GW}^{(2)}$}~ to be consistent with the {\rm NANOGrav}~ result at $2\sigma$ limit.
\begin{figure}[!htbp]
\centering{
\includegraphics[width=0.75\textwidth]{BPLshade_withpbhcont.pdf}
\caption{\small{For a broken power law primordial power spectrum, $A_{\rm CP}$ and $\gamma _{\rm CP}$ fitted for $A_s$ varied in the range $0.002 - 0.02$ for $w=1/3$ (magenta shaded) $w=1$ (cyan shaded), and $w=1/9$ (green shaded). For dark shaded regions, the pivot wavenumber runs from $k_*=1.615\times 10^6$ Mpc$^{-1}$ to $k_*=7.756\times 10^6$ Mpc$^{-1}$ from right to left, where the {\rm NANOGrav}~ fits their observation to a common power law power spectrum. The light shaded regions are shown for the pivot wavenumber $k_*=10^6$ Mpc$^{-1}$ to $k_*=10^7$ Mpc$^{-1}$ from right to left. The dashed contours are for the PBH abundance of $0.1$ (above) and $10^{-5}$ (below) for each case. Black solid contours represent the $1\sigma$ and $2\sigma$ confidence levels for the parameters $A_{\rm CP}$ and $\gamma _{\rm CP}$ observed in the {\rm NANOGrav}~ .}}
\label{BPLcont}}
\end{figure}
}
\section{Discussions and conclusions}
\label{ressec}
The recent {\rm NANOGrav}~ 12.5 yr result can be the first putative signal for the stochastic GW originating in the early universe; further re-analysis of the data will lead to a more precise statement about whether GW has been detected or not. Therefore, it is necessary to check the early universe theories that are consistent with the {\rm NANOGrav}~ result interpreted as GW.

One of the possible explanation for the {\rm NANOGrav}~ result can come from the stochastic GW produced
during inflation. Recently, it was shown in ref.~\cite{Vagnozzi:2020gtf} that if the NANOGrav result is explained by the GW produced during inflation with the standard thermal history, a blue tilted power spectrum will be required. However, the required parameter space to explain the NANOGrav result is entirely excluded by the GW's contribution to effective number of relativistic degrees of freedom during BBN. In this work, we have explored the idea if the {\rm NANOGrav}~ signal can be explained by the GWs generated during inflation with a nonstandard evolution in the post-inflationary cosmological history. While it is well established that BBN must have taken place during the RD era, there is no direct observational probe of EoS for the epoch between the inflation and BBN. In this work, we have assumed a pre-BBN epoch with an equation of state $0\leq w \leq 1$. More specifically, in the first order GW analysis, we have considered three different EoS, $w=0,1/3$, and $1$, which correspond to an early MD era, RD era (standard), and KD era respectively.

\pri{Our main results for the first order GW are summarized in Fig.~\ref{fig:Hinf_nt_allowed_rd_md} where we have shown the {\rm NANOGrav}~ allowed $H_{{\rm inf}}-n_t$ parameter space along with the constraints from CMB, BBN, and LIGO for the three different EoS. It is evident from Fig.~\ref{fig:Hinf_nt_allowed_rd_mdw1} that the parameter space required to explain the {\rm NANOGrav}~ result in a nonstandard thermal history with a KD epoch is ruled out by the bound on $H_{{\rm inf}}$ from the BICEP2/Keck Array and Planck CMB data. It is quite clear from Fig.~\ref{fig:Hinf_nt_allowed_rd_mdw13}  that our results for standard thermal history are consistent with ref~\cite{Vagnozzi:2020gtf} as the entire parameter space required to explain the {\rm NANOGrav}~ result is ruled out by the constraints from BICEP2/Keck Array, Planck CMB, BBN, and LIGO observations. Interestingly, GWs in a nonstandard thermal history with an early MD epoch can explain the {\rm NANOGrav}~ result and can also evade the BICEP2/Keck Array, Planck CMB, BBN and LIGO bounds successfully for a finite region in the parameter space (see Fig.~\ref{fig:Hinf_nt_allowed_rd_mdw0}). This is one of the main results of our analysis for the first order GWs.}

\pri{While calculating the bounds from BBN, the upper limit of the integration given by Eq.~\eqref{eq:bbncons} was set to $10^7$ Hz, which corresponds to a reheat temperature of $10^{15}$ GeV. The upper limit is related to the scale at the end of inflation and lowering its value decreases the value of the integral.} Therefore, for a blue tilted power spectrum, this bound will be further relaxed in case of a lower value of the energy scale at the end of inflation. We have also assumed that the early MD era spans from the end of the inflation until the onset of the RD era at temperature $T_1=10$ MeV. However, there can be a scenario where this MD era lasts for a shorter period, still enveloping the {\rm NANOGrav}~ frequency range. In that case, this bound on GWs from BBN will only be valid if the total contribution of $\Omega_{{\rm GW}}$ from the RD era (lasting in between the end of inflation and the onset of early MD) is smaller than the total contribution from the early MD era.

The constraints from the {\rm NANOGrav}~ on the tilt of the primordial tensor spectrum $n_t$ for different $w$ can also be explained using the slope of the GW spectra. In the first order, the slope of the GW spectrum for a primordial tensor power spectrum with tilt $n_t$ evolving in a $w$ dominated epoch can be calculated from Eq.~\eqref{eq:gammacp} as $5-\gamma _{\rm CP}=\frac{6w-2}{3w+1}+n_t$. Thus, with the $1\sigma$ allowed range $4.5\leq\gamma _{\rm CP}\leq6.5$, the bounds on $n_t$ are: $1-21w\leq 2n_t \leq 5-9w$. Therefore for MD, RD, and KD epochs respectively, the bounds on $n_t$ arising from the {\rm NANOGrav}~ observed slope are: $0.5 \leq n_t \leq 2.5$, $-3 \leq n_t \leq 1$, and $-10 \leq n_t \leq -2$. However, the bounds on the observed amplitude $A_{{\rm CP}}$ at the {\rm NANOGrav}~ along with the bounds from BBN, CMB, and LIGO further reduce these ranges considerably to completely disallow the theoretical signal arising from the first order perturbation theory in the standard RD epoch or an early KD epoch. For an early MD epoch, however, a reduced range of $n_t$ is still consistent with all these bounds as explained in Fig.~\ref{fig:Hinf_nt_allowed_rd_mdw0}.

Our main underlying assumption is that the primordial tensor power spectrum is a pure power law for all the frequencies ranging from CMB scale to the end of inflation. Therefore, in the future, if Laser Interferometer Space Antenna (LISA), which will be observing in the frequency range $f\sim 10^{-2}$ Hz, makes some positive or null detection, that can be used to further constrain the nonstandard pre-BBN epoch. In conclusion, we have explored the possibility of {\rm NANOGrav}~ signal coming from the GWs generated during inflation and evolving in a nonstandard thermal history.  \pri{In the Fig.~\ref{fig:Hinf_nt_allowed_rd_mdw0acp}, we have shown the $\gamma_{{\rm CP}}-A_{{\rm CP}}$ parameter space (blue shaded region) corresponding to the allowed $H_{{\rm inf}}-n_t$ parameter space required to explain the {\rm NANOGrav}~ result in a nonstandard thermal history with an early MD epoch.} If the parameter space  for $\gamma_{{\rm CP}}-A_{{\rm CP}}$ shrinks further with the help of upcoming data from the {\rm NANOGrav}~ or other PTA observatories in the future, it will help in further narrowing down the EoS value of a pre-BBN epoch.

For the second order GW phenomenology in Sec.~\ref{GW2sec}, we consider three separate scenarios for the EoS, standard RD ($w=1/3$), early kinetic energy domination ($w=1$), and a dustlike scenario ($w=1/9$).
\pri{The PBH analysis relevant for such nonstandard scenarios with a general $w$ is discussed in Sec.~\ref{secPBH}. For the primordial power spectrum, we consider three functional forms of $P_{\zeta}(k)$ with (i) a simple Gaussian profile in Sec.~\ref{GPP}, (ii) a lognormal profile with an exponential cut-off in Sec.~\ref{PLN}, and (iii) a broken power law power spectrum in Sec.~\ref{PBPL}. For each of these cases, we quote the typical amplitudes required for each case and then proceed to the relevant implications for PBH.}

For the power spectrum in (i), we find that the analysis to attain a viable scenario that satisfies the {\rm NANOGrav}~ constraints on {$\Omega_{\rm GW}^{(2)}$}~ while also producing abundant PBH crucially depends on not only the amplitude $A_s$ of the primordial curvature power spectrum in Eq.~\eqref{PKG} but also on $\sigma_p$ and $k_*$ for a given $w$. These latter dependences mainly result from the relative shift of both {$\Omega_{\rm GW}^{(2)}$}~ and $\psi_w(M)$ peaks with respect to the peak of $P_{\zeta}(k)$ at $k_*$. For the RD case, both the {$\Omega_{\rm GW}^{(2)}$}~ peak and  $\psi_w(M)$ peak occur very close to $k_*$, therefore the dependence of the PBH abundance and {$\Omega_{\rm GW}^{(2)}$}~ on $\sigma_p$ is mild, as evident from the clustering of the three red points in Fig.~\ref{dotplot}. \pri{For RD, the pivot of the SGWB spectra coincides with the peak scale of the power spectrum i.e., $k_{\rm peak}=k_*$. Since the parameter $A_{\rm CP}$ corresponds to the amplitude of SGWB at the pivot frequency of the NANOGrav, for radiation domination, it only depends on the peak amplitude $A_s$ of the primordial scalar power spectrum. Therefore, the results in the $A_{\rm CP} - \gamma _{\rm CP}$ plane while varying the width $\sigma _p$ (red circle, square, and triangle) are negligibly affected.}

For $w\neq 1/3$, the peak of {$\Omega_{\rm GW}^{(2)}$}~ is far away from $k_*$ and therefore, if the scale of {\rm NANOGrav}~ observation is at $k_0=k_*$ then the peak of {$\Omega_{\rm GW}^{(2)}$}~ is not being observed. Hence, such cases have considerable variations in the $A_{\rm CP}$ - $\gamma _{\rm CP}$ space with change in $\sigma_p$. Now, when we put an additional constraint for obtaining $f_{\rm PBH}=0.1$, we notice the variation for points with different $\sigma_p$ for the two different nonstandard epochs under consideration as the green and blue circles, squares and triangles in Fig.~\ref{dotplot}. This variation shows that given a particular value of $k_*$ in the {\rm NANOGrav}~ observed range, there is a value of $\sigma_p\simeq 2$ where the nonstandard cases come closest to the {\rm NANOGrav}~ $2\sigma$ contour. This feature in the  $A_{\rm CP}$ - $\gamma _{\rm CP}$ space with varying $\sigma_p$ is a result of the competition between the two quantities $\frac{d A_{\rm CP}}{d\sigma_p}$ and $\frac{d f_{\rm PBH}}{d\sigma_p}$. It is analytically challenging to check the variation of these two terms since both of them involve rigorous numerical integrations.

Now, fixing $\sigma_p=2$ and varying $A_s$ as per Table~\ref{TT1} and $k_*$ in the {\rm NANOGrav}~ range, we show the contours (dashed lines in Fig.~\ref{regcont}) for abundant PBH production or less for the power spectrum (i). We conclude from Fig.~\ref{regcont} that the $w=1$ case cannot produce a copious amount of PBH and satisfy the {\rm NANOGrav}~ constraints simultaneously. We notice that for the $w=1/9$ case, the contour for $f_{\rm PBH}=0.1$ just enters the {\rm NANOGrav}~ $2\sigma$ contour in this range of $k_*$. Guided by our previous understanding of the relative shift in peak positions, we find that for $w=1/9$, choosing a different $k_*$ range such that {\rm NANOGrav}~ observes closer to the peak of the {$\Omega_{\rm GW}^{(2)}$}~ spectrum makes the nonstandard $w=1/9$ case consistent with the {\rm NANOGrav}~ $2\sigma$ contour. The gain in a $w=1/9$ case over a pure RD epoch is that both abundant PBH and {\rm NANOGrav}~ consistency are reached with a smaller $A_s\sim 0.007$ in the $w=1/9$ case as compared to $A_s\sim 0.01$ required in the RD case.

The reason for the consistency of the $w=1/9$ case and the inconsistency of the $w=1$ case for the Gaussian power spectrum can also be explained in terms of the slope $\gamma _{\rm CP}$. In a general $w$-dominated epoch, the slope of the {$\Omega_{\rm GW}^{(2)}$}~ spectrum is $5-\gamma _{\rm CP}=\frac{6w-2}{3w+1}+2n_2$, where $n_2$ is the slope of $P_{\zeta}(k)$. For a Gaussian primordial power spectrum given in Eq.~\eqref{PKG}, $n_2=-\frac{1}{\sigma_p^2}\log_{10}(k/k_*)$. Taking the $1\sigma$ limit on $\gamma _{\rm CP}$ from the {\rm NANOGrav}~ $4.5\lesssim \gamma _{\rm CP}\lesssim 6.5$, we find no real limit for $\sigma_p $ for the $w=1$ case. However, for the $w=1/9$ case, the lower limit on $\gamma _{\rm CP}$ gives $\sigma_p \gtrsim 1.155$ with $k=k_0=3.6\times 10^6$ Mpc$^{-1}$ and $k_*\simeq 3.6\times 10^7$ Mpc$^{-1}$. Hence, the scenario where a nonstandard $w=1/9$ epoch is consistent with {\rm NANOGrav}~ result and abundant PBH production (green dashed contour in Fig.~\ref{regcont2}), $n_2 \vert_{k_0}\lesssim 0.5$.

\pri{The scalar power spectrum for the case (ii) is necessarily a tilted Gaussian in the $k$-space, and therefore, we expect our results for this case to be similar to the symmetric Gaussian power spectrum in case (i). Form Fig.~\ref{LNcont}, it is clear that $w=1$ case results in insignificant PBH abundance for the predicted {$\Omega_{\rm GW}^{(2)}$}~ to be consistent with the {\rm NANOGrav}~ data; the $w=1/3$ scenario is consistent with the {\rm NANOGrav}~ data within $2\sigma$ confidence limit for $f_{\rm PBH}\simeq 0.1$; the $w=1/9$ case is consistent with the {\rm NANOGrav}~ data for $f_{\rm PBH}> 0.1$ only for $k_* > 5\times 10^6$ Mpc$^{-1}$. Indeed, all of these results are similar to the results for the symmetric Gaussian case in Sec.~\ref{GPP}. The difference between the exact results for {$\Omega_{\rm GW}^{(2)}$}~ and PBH abundance between Sec.~\ref{GPP} and Sec.~\ref{PLN} is mainly due to the difference in their width.
For the $P_{\zeta}(k)$ considered in Sec.~\ref{PLN}, the advantage for the $w=1/9$ epoch over typical RD epoch is that the {\rm NANOGrav}~ consistency and $f_{\rm PBH}\gtrsim 0.1$ can be simultaneously satisfied for $A_s(w=1/9)=0.0078$ as compared to that required for the RD, i.e., $A_s(w=1/3)=0.011$.}

\pri{The third power spectrum considered in Sec.~\ref{PBPL} is a broken power law, for which the implications for PBH and second order GW are given in Fig.~\ref{BPLcont} in a consolidated manner. In this case, out of the three post-inflationary evolutions considered, abundant PBH formation with $f_{\rm PBH}=0.1$ is consistent with the {\rm NANOGrav}~ result only for RD. For the nonstandard cases, {$\Omega_{\rm GW}^{(2)}$}~ in the {\rm NANOGrav}~ observed range leads to a tiny value of $f_{\rm PBH}$ (for $w=1$) or it overshoots the theoretical limit $f_{\rm PBH}=1$ (for $w=1/9$). The difference in conclusion for the  $w=1/9$ case for this power spectrum as compared to the other two $P_{\zeta}(k)$ considered in Sec.~\ref{GPP} and Sec.~\ref{PLN} can be attributed to the sharp drop in power for the green curve in Fig.~\ref{OmgA2_BPL} from the peak to $k=k_*$.}

\begin{figure}[!htbp]
\centering{
\includegraphics[width=0.75\textwidth]{OMG2fullplot.pdf}
\caption{\small{The second order GW spectra for the three values of $w$ for the Gaussian power spectrum in Eq.~\ref{PKG} (solid lines), lognormal power spectrum with an exponential cutoff in Eq.~\ref{PKLN} (dashed lines), and broken power law power spectrum in Eq.~\ref{PKBPL} (dotted lines) are shown. The red, blue and green curves represent $w=1/3$, $1$, and $1/9$ respectively. LIGO runs for O1, O2, and O5 (future) are given as cyan curves (top to bottom); magenta shade represnents the proposed sensitivity for the LISA survey; black shade region represent $1\sigma$ bound from the {\rm NANOGrav}~ 12.5 year observation; orange shade region represents the proposed sensitivity from the future SKA survey.}}
\label{fullomg2obs}}
\end{figure}
\pri{Similar to Table~\ref{TT2}, PBH mass ranges with different choices of $k_*$ for the second and third form of $P_{\zeta}(k)$ and for all three $w$ under consideration are presented in Table~\ref{TT4}.}
\begin{table}[!htbp]
\caption{Ranges of PBH mass $M$ produced in abundance for each EoS}
\begin{center}
\label{TT4}
\resizebox{0.85\textwidth}{!}{
\begin{tabular}{|c|c|c|c|c|}
\hline
$P_{\zeta}(k)$ & $k_*$ in Mpc$^{-1}$ & $M/M_{\odot} (w=1/3)$ & $M/M_{\odot} (w=1)$ & $M/M_{\odot} (w=1/9)$ \\
\hline
\hline
Lognormal with exponential cutoff & $10^6$  & $0.8-40$ & $2-200$ &  $0.08 - 100$\\
\hline
& $10^7$  & $7\times 10^{-3} - 0.5$ & $0.07 - 8$ & $5\times 10^{-4} - 0.5$\\
\hline
Broken power law & $10^6$  & $0.2 - 10$ & $1 - 30$ & $0.02 - 8$\\
\hline
& $10^7$  & $2\times 10^{-3} - 0.2$ & $0.025 - 1$  &  $10^{-4} - 0.03$\\
\hline
\end{tabular}}
\end{center}
\end{table}

The discussions in terms of the slopes of the second and first order GW spectrum in this section emphasise its significance in constraining pre-BBN cosmologies and primordial power spectra when a (probable) GW detection is available. We have checked that our results satisfy BBN bound on {$\Omega_{\rm GW}^{(2)}$}~ given in Eq.~\eqref{eq:bbncons}. The range of possible PBH masses is different for different $w$ and $k_*$ values and different $P_{\zeta}(k)$ (see Tables~\ref{TT2} and~\ref{TT4}). In case of the Gaussian power spectrum, for the nonstandard case $w=1/9$ that is finally consistent in our analysis can generate PBH in the mass range $5\times 10^{-6}M_{\odot}$ to $3\times 10^{-3}M_{\odot}$ depending on $k_*$. The observational constraints from lensing experiments such as Subaru HSC~\cite{Niikura:2017zjd}, EROS/MACHO~\cite{Tisserand:2006zx} and OGLE~\cite{Niikura:2019kqi} put upper bound on $f_{\rm PBH}$ to be 0.1 in a few points across this mass range. However, the $f_{\rm PBH}\simeq 10^{-3}$ contour will still be inside {\rm NANOGrav}~ $2\sigma$ in this case, which is observationally valid throughout this mass range. For the RD epoch, near solar mass PBHs can be generated with $f_{\rm PBH}=0.1$, which is the maximum observationally allowed abundance in this range, constrained mainly with EROS/MACHO~\cite{Tisserand:2006zx}.

\pri{For the lognormal power spectrum in Eq.~\ref{PKLN}, the minimum and maximum masses of PBH that can be produced for the standard RD case are $7\times 10^{-3}M_{\odot}$ and $40M_{\odot}$ respectively, which can be compared with the observational bounds discussed above in a similar manner. For this $P_{\zeta}(k)$ with $w=1/9$, both abundant PBH formation and {\rm NANOGrav}~ consistency are satisfied only for $k_*=10^7$ Mpc$^-1$ in Table~\ref{TT4}, where PBH masses are produced in the range $5\times 10^{-4} - 0.5 M_{\odot}$. Observational constraints on PBH abundance from EROS, OGLE, Icarus and Microlensing constaints from supernova are effective in this mass range and maximum possible abundance is $f_{\rm PBH}< 0.1$. For the broken power law power spectrum in Eq.~\ref{PKBPL}, the only viable case is RD, which can lead to abundant PBH production in the mass range $2\times 10^{-3}- 10 M_{\odot}$ and is constrained by the experiments discussed above constraining $f_{\rm PBH}\simeq 0.1$. Therefore, with the choice of parameters and $k_*$, for all of the $P_{\zeta}(k)$ under consideration, $f_{\rm PBH}= 0.1$ is a good benchmark to explore, which is shown as dashed contours in Fig.s~\ref{regcont}, ~\ref{regcont2}, ~\ref{LNcont} and~\ref{BPLcont}. As pointed out in earlier sections, the mass ranges for all theree $P_{\zeta}(k)$ discussed here only envelop a fraction of the mass range where LIGO/Virgo survey has found binary merger events. For the cases where {\rm NANOGrav}~ observation is consistent with abundant PBH formation (see Table~\ref{TT3}) can only lead to very small PBH abundance in the LIGO/Virgo mass range. Therefore, in this set up, it is unlikely that these the BH mergers observed in LIGO/Virgo surveys are of primordial nature~\cite{Bird:2016dcv,Clesse:2016vqa,Sasaki:2016jop}.}

\pri{The full {$\Omega_{\rm GW}^{(2)}$}~$h^2$ for the three forms of $P_{\zeta}(k)$ and three values of $w$ are plotted in Fig.~\ref{fullomg2obs} with a few present and future observational bounds. In this plot, $A_s =0.007$ for all the cases, $\sigma _p=2$ for the Gaussian power spectrum (solid lines); $\alpha=0.17$ and $\beta=0.62$ for the lognormal power spectrum with an exponential cutoff (dashed lines); and $\alpha=4$ and $\beta=0.6$ for the broken power law power spectrum in (dotted lines). The red, blue and green curves represent $w=1/3$, $1$, and $1/9$ respectively. Evidently, future Square Kilometre Array (SKA)~\cite{Zhao:2013bba} data will be crucial in constraining the {$\Omega_{\rm GW}^{(2)}$}~ spectrum in this frequency range. The conclusions for different $w$ and different $P_{\zeta}(k)$ for $k_*=3.6\times 10^6$ Mpc$^{-1}$ are summarised in Table~\ref{TT3}. Evidently, modifying the values of the parameters $\alpha$ and $\beta$ for power spectra in Sec.~\ref{PLN} and~\ref{PBPL} will modify the conclusions of this analysis. A straightforward way to check the dependence of this analysis on these parameters is to perform a similar analysis as done in Fig.~\ref{regcont2} of Sec.~\ref{GPP}, which we hope to explore in future studies.}

\begin{table}[!htbp]
\caption{Summary: consistency of PBH formation with $f_{\rm PBH}=0.1$ while satisfying {\rm NANOGrav}~ 12.5 year data as second order GW spectra for different $P_{\zeta}(k)$. }
\begin{center}
\label{TT3}
\resizebox{0.95\textwidth}{!}{
\begin{tabular}{|c|c|c|c|c|}
\hline
$P_{\zeta}(k)$ & $k_*$ (Mpc$^{-1}$) & $w=1/3$ & $w=1$ & $w=1/9$ \\
\hline
\hline
Gaussian  & $10^6$ & consistent at $2\sigma$ with $A_s\simeq 0.01$ & inconsistent & inconsistent ($f_{\rm PBH}>1$)\\
\hline
 & $10^7$ & consistent at $2\sigma$ with $A_s\simeq 0.01$ & inconsistent & consistent with $A_s\simeq 0.007$\\
\hline
Lognormal with exponential cutoff & $10^6$ & inconsistent at $2\sigma$ with $A_s\simeq 0.011$ & inconsistent & inconsistent (large $\gamma _{\rm CP}$) \\
\hline
  & $10^7$ & consistent at $2\sigma$ with $A_s\simeq 0.011$ & inconsistent & consistent with $A_s\simeq 0.0078$ \\
\hline
Broken power law & $10^6$ & consistent at $1\sigma$ with $A_s\simeq 0.012$ & inconsistent & inconsistent ($f_{\rm PBH}>1$) \\
\hline
 & $10^7$ & consistent at $2\sigma$ with $A_s\simeq 0.0122$ & inconsistent & inconsistent ($f_{\rm PBH}>1$)\\
\hline
\end{tabular}}
\end{center}
\end{table}

After the announcement of the {\rm NANOGrav}~ 12.5 yr data release, many works have studied the consistency of abundant PBH formation with the {\rm NANOGrav}~ observation in terms of {$\Omega_{\rm GW}^{(2)}$}~ . References~~\cite{Vaskonen:2020lbd,DeLuca:2020agl,Kohri:2020qqd} have studied this in RD epochs. Ref.~\cite{DeLuca:2020agl} considered a wide constant primordial power spectrum and found high abundance in the PBH mass range $10^{18}-10^{21}$ gm, while still satisfying the {\rm NANOGrav}~ constraints. Ref.~\cite{Vaskonen:2020lbd} studied different primordial curvature power spectra. The recent work~\cite{Domenech:2020ers} discusses implications of the {\rm NANOGrav}~ result on the formation of planetary mass PBHs in nonstandard dust-like epochs with $-0.091<w<0.048$. They use narrow peak for their analysis, specifically a Dirac delta function, and translate the constraints on {$\Omega_{\rm GW}^{(2)}$}~ to that of the amplitude which generates a monochromatic mass function for PBH, whereas our $P_{\zeta}(k)$ with $\sigma_p=2$ can be noted as a broad spectrum.

\pri{Ref.~\cite{Kohri:2020qqd} motivates a narrow Gaussian curvature power spectrum using the Press-Schechter formalism to generate abundant near solar mass PBHs to be consistent with the {\rm NANOGrav}~ results in RD. In our work with the Gaussian power spectrum in Sec.~\ref{GPP}, we also consider Press-Schechter formalism for PBH analysis and found out $\sigma_p=1$ case to be consistent with~\cite{Kohri:2020qqd} for RD.
The lognormal power spectrum with an exponential cut off in Sec.~\ref{PLN} and the broken power law power spectrum in Sec.~\ref{PBPL} has also been studied in~\cite{Vaskonen:2020lbd} for RD epoch, however, they had discarded these scenarios because its {\rm NANOGrav}~ consistent amplitude does not work to generate abundant PBH. However, in our case, abundant PBH production in RD for both of these power spectra are consistent with the {\rm NANOGrav}~ data and therefore disagrees with the results presented in~\cite{Vaskonen:2020lbd}. The reasons for this disagreement are similar to those pointed out in the discussion section of~~\cite{Kohri:2020qqd}. Whereas ref.~\cite{Vaskonen:2020lbd} consider critical collapse formalism to evaluate PBH abundance, we note that it is only valid for exact spherical collapse, and we do not consider this formalism for PBH analysis in this paper. Other reasons for the disagreement include a choice of wider window functions for PBH analysis in our case; inclusion of the transfer function in the definition of $\sigma ^2 _\delta$, and inclusion of the nonlinear relation between density perturbation and curvature perturbation in ref.~\cite{Vaskonen:2020lbd}.}

For a broad primordial power spectrum in general, the resulting {$\Omega_{\rm GW}^{(2)}$}~ spectrum will also be broad. Therefore, the future GW detectors close to the {\rm NANOGrav}~ frequency range will further constrain the width through (non-)observation of the stochastic GW background. These future constraints combined with the {\rm NANOGrav}~ limits will motivate a thorough understanding of the primordial (inflationary) and pre-BBN cosmologies and constrain the PBH abundances more stringently.

For our case the mass function is extended, including both solar mass and planetary mass PBHs. The PBH analysis also depends on the formalism under consideration so that abundances from peak theory analysis or critical collapse analysis may differ considerably from the Press-Schechter analysis used here. We have used a broad window function for the Gaussian power spectrum, which substantially helps in enhancing PBH. We note here that for a broad power spectrum, the integration kernels have to be evaluated on a case by case basis for different EoS and therefore, it is challenging to show constraints as a functions of $w$ directly. In a more explicit scenario, it is possible to start from a well-motivated theoretical framework where dust-like nonstandard epochs can arise and deriving a thorough evolution for $\Omega_{\rm GW}^{(1)}$ and {$\Omega_{\rm GW}^{(2)}$}~ in this nonstandard epoch and RD, it is possible to constrain parameters of the theoretical interest using the {\rm NANOGrav}~ result and upcoming observations of stochastic GW background. Also, the {$\Omega_{\rm GW}^{(2)}$}~ for pure matter dominated case $w=0$ should be calculated separately for a broad primordial power spectrum since one needs to take care of the scales where non-linear growths become important. For $w=0$, the PBH formation mechanism is also quite different from $w>0$ cases due to the absence of pressure and often results in PBHs with nonzero-spins and aspherical shapes. We hope to come back to these issues in the near future.

\section*{Acknowledgement}
The authors sincerely thank the referees for their insightful comments and suggestions that helped in improving the structure and quality of this paper. S.B. is supported by the institute postdoctoral fellowship from Physical Research Laboratory, India.

\bibliographystyle{JHEP}