# The LHCb state $P_{\psi s}^{\Lambda}(4338)$ as a triangle singularity

T. J. Burns  Department of Physics, Swansea University, Singleton Park,
Swansea, SA2 8PP, UK.  E. S. Swanson  Department of Physics and Astronomy,
University of Pittsburgh, Pittsburgh, PA 15260, USA.

###### Abstract

We present a model for the $J/\psi\,\Lambda$ spectrum in $B^{-}\to
J/\psi\,\Lambda\,\bar{p}$ decays, including the $P_{\psi s}^{\Lambda}(4338)$
baryon recently observed by the LHCb collaboration. We assume production via
triangle diagrams which couple to the final state via non-perturbative
interactions which are constrained by heavy-quark and $SU_{3}$-flavor
symmetry. The bulk of the distribution is described by a triangle diagram with
a color-favored electroweak vertex, while the sharp $P_{\psi
s}^{\Lambda}(4338)$ enhancement is due to the triangle singularity in another
diagram featuring a $1/2^{-}$ baryon consistent with $\Sigma_{c}(2800)$. We
predict a comparable $P_{\psi s}^{\Lambda}(4338)$ signal in
$\eta_{c}\,\Lambda$, and anticipate possible large isospin mixing effects
through decays to $J/\psi\,\Sigma^{0}$ and $\eta_{c}\,\Sigma^{0}$.

##  I Introduction

The LHCb collaboration continues its exploration of hadronic interactions as
revealed by electroweak decays of heavy hadrons, recently announcing the
discovery of a signal in the $J/\psi\,\Lambda$ mass spectrum of $B^{-}\to
J/\psi\,\Lambda\,\bar{p}$ ChenParticleZoo . The mass and width of $P_{\psi
s}^{\Lambda}(4338)$ are

$\displaystyle M$ $\displaystyle=4338.2\pm 0.7~{}\textrm{ MeV},$ (1)
$\displaystyle\Gamma$ $\displaystyle=7.0\pm 1.2~{}\textrm{ MeV},$ (2)

and $J^{P}=1/2^{-}$ quantum numbers are preferred.

Because of its proximity to $\Xi_{c}\bar{D}$ threshold, a molecular
interpretation of $P_{\psi s}^{\Lambda}(4338)$ has been proposed
Karliner:2022erb ; Wang:2022mxy ; Yan:2022wuz . Molecules with
$\Xi_{c}\bar{D}$ constituents have been predicted in a wide range of models,
typically assuming a binding interaction due to boson exchange, or effective
field theory constrained by heavy quark symmetry Wang:2019nvm ; Xiao:2019gjd ;
Peng:2019wys ; Peng:2020hql ; Chen:2020uif ; Dong:2021juy ; Yan:2021nio ;
Chen:2021cfl ; Chen:2020kco ; Liu:2020hcv ; Zhu:2021lhd ; Xiao:2021rgp ;
Chen:2022onm ; Wang:2019nvm . In such models, whether or a not a particular
state binds is ultimately determined by one or more parameters which are fit
to data, such as the form factor cut-off, or contact terms attributed to
unknown short-distance physics. For this reason, the molecular approach is
only robust to the extent that it can simultaneously describe several
different states with the same fit parameters. With respect to this criterion,
the molecular scenario for $P_{\psi s}^{\Lambda}(4338)$ runs into problems.

If $P_{\psi s}^{\Lambda}(4338)$ is a $\Xi_{c}\bar{D}$ molecule ($1/2^{-}$), it
could potentially have $\Xi_{c}\bar{D}^{*}$ partners ($1/2^{-}$ and
$3/2^{-}$). From heavy quark symmetry, the potentials in all three channels
are identical,

$\displaystyle
V(\Xi_{c}\bar{D},1/2^{-})=V(\Xi_{c}\bar{D}^{*},1/2^{-})=V(\Xi_{c}\bar{D}^{*},3/2^{-}),$
(3)

which implies the three states should have comparable binding energies, an
expectation which is confirmed in a one-boson model calculation Wang:2022mxy .
The $P_{cs}(4459)$ state observed at LHCb Aaij:2020gdg  is a candidate for a
$\Xi_{c}\bar{D}^{*}$ molecule Peng:2020hql ; Chen:2020kco ; Liu:2020hcv ;
Zhu:2021lhd ; Xiao:2021rgp , but if it is the partner of $P_{\psi
s}^{\Lambda}(4338)$ as a $\Xi_{c}\bar{D}$ molecule, it implies a drastic
violation of heavy-quark symmetry, since its binding energy is so large (19
MeV) compared to $P_{\psi s}^{\Lambda}(4338)$ (which is at threshold).

Similarly, the hypothesis that $P_{cs}(4459)$ consists of two states
Aaij:2020gdg , interpreted as $\Xi_{c}\bar{D}^{*}$ molecules ($1/2^{-}$ and
$3/2^{-}$) Karliner:2022erb ; Wang:2022mxy , is problematic, because the 13
MeV mass splitting contradicts the above expectation from heavy-quark
symmetry. Although some level of splitting is expected in models (for example
due to coupled-channel effects) Wang:2019nvm ; Peng:2020hql , it remains to be
seen whether such models can obtain the right level of splitting, and also
reconcile the significant $\Xi_{c}\bar{D}^{*}$ binding with the lack of
binding in $\Xi_{c}\bar{D}$. We note that in the model of ref. Yan:2022wuz ,
there is no region of parameter space in which both $P_{\psi
s}^{\Lambda}(4338)$ and two $P_{cs}(4459)$ states can be accommodated.

Refs. Karliner:2022erb ; Wang:2022mxy  have also argued for an analogy between
$P_{\psi s}^{\Lambda}(4338)$ and $P_{cs}(4459)$ (as $\Xi_{c}\bar{D}^{(*)}$
molecules), and $P_{c}(4312)$, $P_{c}(4440)$ and $P_{c}(4457)$ Aaij:2015tga ;
Aaij:2019vzc ; Aaij:2016phn  (as $\Sigma_{c}\bar{D}^{(*)}$ molecules). The
analogy is misleading, considering that the corresponding potentials are
neither related by $SU(3)$ flavor (as $\Xi_{c}$ and $\Sigma_{c}$ belong to
different flavor multiplets), nor heavy-quark spin symmetry (as $\Xi_{c}$ and
$\Sigma_{c}$ have different light quark spins) Peng:2020hql . Indeed, heavy
quark symmetry implies a completely different pattern of states in
$\Sigma_{c}\bar{D}^{(*)}$ systems (where the potentials are spin-dependent
Liu:2019tjn ; Liu:2019zvb ; Du:2019pij ; Sakai:2019qph ; Valderrama:2019chc ;
Xiao:2019aya ; Peng:2020xrf ; Du:2021fmf ; Liu:2020hcv ) and
$\Xi_{c}\bar{D}^{(*)}$ systems (where they are not). Moreover, the analogy
relies on the assumption that $P_{c}(4440)$ and $P_{c}(4457)$ are both
$\Sigma_{c}\bar{D}^{*}$ molecules, and we recently argued that this assumption
is not consistent with experimental constraints Burns:2021jlu . (Scenarios
with different interpretations for $P_{c}(4457)$ Burns:2019iih ; Burns:2022uiv
do not have the same problem.)

An additional awkward feature of the molecular scenario for $P_{\psi
s}^{\Lambda}(4338)$, which has hardly been discussed in the literature, is
that it is not bound with respect to the $\Xi_{c}\bar{D}$ thresholds, but
somewhat above:

$\displaystyle\Xi_{c}^{0}\bar{D}^{0}$ $\displaystyle=4335.28\pm
0.33~{}\textrm{ MeV},$ (4) $\displaystyle\Xi_{c}^{+}D^{-}$
$\displaystyle=4337.37\pm 0.28~{}\textrm{ MeV}.$ (5)

The situation here is similar to $P_{c}(4457)$, which is widely interpreted as
a $\Sigma_{c}\bar{D}^{*}$ bound state, despite having a mass which is
consistent with the threshold not only for $\Sigma_{c}\bar{D}^{*}$, but also
$\Lambda_{c}(2595)\bar{D}$. Signals at (rather than below) threshold are more
amenable to non-resonant interpretations, and we recently demonstrated that
$P_{c}(4457)$ can be explained as a cusp or an enhancement due to the
logarithmic triangle singularity Burns:2022uiv . In this paper we explore
related possibilities for $P_{\psi s}^{\Lambda}(4338)$.

##  II Model

We assume, as in our previous work Burns:2020epm ; Burns:2022uiv , that the
distribution can be described by triangle diagrams which couple to the final
state through interactions which respect heavy-quark symmetry, and that the
dominant diagrams are those with color-favored weak vertices. Hence we
consider $\bar{B}\to D_{s}^{(*)-}\bar{D}^{(*)}$ transitions, noting that such
branching fractions range from approximately 1-3%. We can then form the
triangle diagram shown in Fig. 1 (left), involving virtual $\Lambda_{c}$
exchange, followed by rescattering into the final $J/\psi\,\Lambda$.
Generically, the distribution associated with this diagram peaks around the
$\Lambda_{c}^{+}D_{s}^{-}$ threshold (where it has a cusp). We notice that the
$J/\psi\Lambda$ distribution ChenParticleZoo  has exactly this shape, and we
regard this as strong support for this proposed production mechanism. We also
notice that, similar to the other exotic hadron systems Burns:2021jlu ;
Burns:2022uiv , the tree-level diagrams for the $J/\psi\,\Lambda\,\bar{p}$
final state are color-suppressed; hence it is natural to assume that the
color-favored triangle diagram is a dominant contribution.

Figure 1: Diagram $T_{1}$ (left) has a color-favored weak transition, whereas
$T_{2}$ (right) is color-suppressed but enhanced at the $\Xi_{c}\bar{D}$
threshold due to the triangle singularity. The solid circles indicate non-
perturbative final-state interactions, as described in the text.

The proximity of $P_{\psi s}^{\Lambda}(4338)$ to $\Xi_{c}\bar{D}$ threshold
also suggests a role for diagrams in which $\Xi_{c}\bar{D}$ rescatters to
$J/\psi\,\Lambda$, and here we notice an intriguing possibility. Such a novel
intermediate state can be realized via an electroweak decay, such as
$B^{-}\to\bar{\Lambda}_{c}\Xi_{c}$ or $B^{-}\to\bar{\Sigma}_{c}\Xi_{c}$, as
shown in Fig. 1 (right). Although the electroweak vertex is color-suppressed,
we notice that the first of these modes has been observed in experiment
Belle:2018kzz , with quite significant branching fraction $(9.51\pm 2.1\pm
0.88)\times 10^{-4}$. Even accounting for some suppression at the electroweak
vertex, such diagrams can make a significant contribution to the amplitude
near the $\Xi\bar{D}$ threshold, if the mass of the $\bar{\Lambda}_{c}$ or
$\bar{\Sigma}_{c}$ leads to a logarithmic singularity in the triangle diagram.
The same idea has been applied in different contexts to explain other exotic
hadron phenomena, as reviewed in ref. Guo:2019twa .

We will concentrate on the $\bar{\Sigma}_{c}$ (rather than
$\bar{\Lambda}_{c}$) diagram, since only this is capable of producing the
specific charge combination $\Xi_{c}^{+}D^{-}$ which, given the mass of
$P_{\psi s}^{\Lambda}(4338)$, appears to be relevant. Solving the non-
relativistic version of the Landau equations Guo:2019twa , we find that the
requisite $\bar{\Sigma}_{c}$ mass is 2810 MeV. Amazingly, there is a three-
star resonance $\Sigma_{c}(2800)$ whose mass is almost exactly right. The
skeptical reader may suspect that there is a large number of possible
$\Sigma_{c}$ states, and that we are simply choosing one with the right mass
in order to make our proposed mechanism work. On the contrary – aside from the
familiar ground states $\Sigma_{c}(2455)$ and $\Sigma_{c}(2520)$, the only
$\Sigma_{c}$ baryon in the Particle Data Group tables is $\Sigma_{c}(2800)$
ParticleDataGroup:2020ssz .

Hence we will attempt to fit the $J/\psi\,\Lambda$ distribution, adopting the
following model for the $B^{-}\to J/\psi\,\Lambda\,\bar{p}$ amplitude:

$\mathcal{A}=b+g_{1}T_{1}+g_{2}\frac{1}{\sqrt{6}}\left[2T_{2}^{(--)}-T_{2}^{(-)}\right],$
(6)

where $b$ (the background, a complex constant) and $g_{1,2}$ (the production
couplings) are fit to data, and $T_{1}$ and $T_{2}$ are the sub-amplitudes
corresponding to the diagrams in the left and right panels of Fig. 1,
respectively, computed as described below.

Because the $P_{\psi s}^{\Lambda}(4338)$ peak is much closer to the threshold
for $\Xi_{c}^{+}D^{-}$ rather than $\Xi_{c}^{0}\bar{D}^{0}$, we do not expect
isospin symmetry to be respected in this system. This was discussed in refs
Yan:2022wuz ; Meng:2022wgl , and is similar to related observations in, for
example, $X(3872)$ and the $P_{c}$ states Swanson:2003tb ; Close:2003sg ;
Burns:2015dwa ; Guo:2019fdo ; Guo:2019kdc . In $T_{2}$ we therefore consider
the $\Xi_{c}^{+}D^{-}$ and $\Xi_{c}^{0}\bar{D}^{0}$ diagrams separately, and
refer to these as $T_{2}^{(--)}$ and $T_{2}^{(-)}$, corresponding to
$\bar{\Sigma}_{c}^{--}$ and $\bar{\Sigma}_{c}^{-}$ states in the triangle
diagrams. The remaining factors in equation (6) arise from an isospin
decomposition of the amplitude. Notice that the weighted combination of
$T_{2}^{(--)}$ and $T_{2}^{(-)}$ corresponds to $\Xi_{c}\bar{D}$ in a linear
combination of isospin 0 and 1. Hence the production mechanism itself does not
respect isospin, even before considering the mass difference between
$\Xi_{c}^{+}D^{-}$ and $\Xi_{c}^{0}\bar{D}^{0}$. This is a different (and more
significant) source of isospin mixing than that due to the $\Xi_{c}\bar{D}$
masses; similar effects could be present in the $X(2900)$ system Burns:2020xne
.

Our calculation of the amplitudes $T_{1}$ and $T_{2}$ follows the method
outlined in detail in our previous paper Burns:2022uiv , so here we just
mention some key points. The triangle diagrams are computed in the
nonrelativistic limit, incorporating form factors for the strong decay vertex
and the final state interactions. The form factor scale is fixed at 800 MeV,
as in the previous study. Because modeling the electroweak vertex is
difficult, and because it does not vary much over the rather narrow phase
space available, we have chosen to employ a constant electroweak vertex, whose
strength is ultimately absorbed into production couplings $g_{1}$ and $g_{2}$
which, as discussed below, are fit to data.

In diagram $T_{1}$, the functional dependence of the amplitude is insensitive
to the choice of the “$\bar{D}^{(*)}$” mass. Hence instead of separately
considering contributions from $\bar{D}$ and $\bar{D}^{*}$, we compute a
single diagram, using the physical $\bar{D}^{*}$ mass.

In diagram $T_{2}$, we notice that, in order to generate a prominent
enhancement at $\Xi_{c}\bar{D}$ threshold, we need the
$\bar{\Sigma}_{c}\to\bar{D}\bar{p}$ vertex to be S-wave. The corresponding
$\Sigma_{c}$ state therefore has $J^{P}$ quantum numbers $1/2^{-}$. The
quantum numbers of $\Sigma_{c}(2800)$ have not been measured, but the mass is
consistent with expectations for the 1P multiplet, which includes states with
quantum numbers $1/2^{-}$, $3/2^{-}$ and $5/2^{-}$ Copley:1979wj ; Chen:2007xf
; Cheng:2006dk ; Cheng:2015naa ; Wang:2021bmz . We therefore assume that there
is a $1/2^{-}$ state around 2800 MeV, though we are not necessarily assuming
it is $\Sigma_{c}(2800)$ itself. For the purposes of the calculation, we fix
the “$\bar{\Sigma}_{c}$” mass to 2801 MeV (the measured mass of
$\Sigma_{c}(2800)^{++}$ ParticleDataGroup:2020ssz ), but we notice that the
fit quality is not highly sensitive to this choice. (Although the logarithmic
singularity strictly arises within a specific and narrow window of
$\bar{\Sigma}_{c}$ masses, in practice we find very strong enhancements in the
triangle diagram over a range of $\bar{\Sigma}_{c}$ masses around 2800 MeV.)

Most particles in the triangle diagrams are very narrow, so we ignore their
widths and instead introduce a small imaginary part $\epsilon$ in the energy
denominators. The exception is $\bar{\Sigma}_{c}$, whose width is included
explicitly. We consider two cases: $\Gamma=70$ MeV (from the measured width of
the $\Sigma_{c}(2800)$), and $\Gamma=15$ MeV (a model prediction for the width
of the $1/2^{-}$ state in the 1P multiplet Wang:2021bmz ).

The triangle diagrams couple to the $J/\psi\,\Lambda$ final state via non-
perturbative final-state interactions, represented by the solid circles in
Fig. 1. Following ref. Burns:2022uiv , we assume a separable form for the
interaction potential, and so obtain the non-perturbative T-matrix by solving
the Bethe-Heitler equation, $T=V+VGT$, using algebraic methods. The final-
state interactions are responsible for couplings among all the relevant
channels so far discussed ($\Lambda_{c}^{+}D_{s}^{-}$, $\Xi_{c}^{+}D^{-}$,
$\Xi_{c}^{0}\bar{D}^{0}$, $\Lambda J/\psi$), as well as others which are
related to these by heavy-quark symmetry and which also couple to $1/2^{-}$ in
S-wave. In particular, we include $\Lambda\eta_{c}$, $\Sigma J/\psi$ and
$\Sigma\eta_{c}$ as possible final states of interest, noting that the
$\Sigma$ modes are relevant because of the explicit isospin mixing in the
model. We do not include other channels such as $\Lambda_{c}D_{s}^{*}$,
$\Xi_{c}\bar{D}^{*}$ and $\Xi_{c}^{(\prime,*)}\bar{D}^{(*)}$, whose thresholds
are beyond the kinematic boundary for $J/\psi\Lambda$ in $B^{-}\to
J/\psi\,\Lambda\,\bar{p}$.

We choose to model the final-state interactions as contact terms constrained
by heavy-quark symmetry, following an approach which has been widely applied
to $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ systems Liu:2018zzu ; Liu:2019tjn ;
Sakai:2019qph ; Valderrama:2019chc ; Liu:2019zvb ; Du:2019pij ; Peng:2020xrf ;
Du:2021fmf , and more recently $\Xi_{c}^{(\prime,*)}\bar{D}^{(*)}$ systems
Peng:2019wys ; Yan:2021nio ; Peng:2020hql ; Liu:2020hcv ; Yan:2022wuz . We
previously tabulated the relevant contact terms for S-wave interactions among
isospin 1/2 channels $\Lambda_{c}\bar{D}^{(*)}$,
$\Sigma_{c}^{(*)}\bar{D}^{(*)}$, $NJ/\psi$ and $N\eta_{c}$ Burns:2022uiv .
Most of the contact terms we need for the present case can be extracted from
those by assuming, as in refs Peng:2020hql ; Yan:2021nio ; Yan:2022wuz , that
the interactions are invariant under rotations in $SU_{3}$ flavor space. Hence
the matrix elements for the
$\big{|}SU_{3}\textrm{~{}flavor},SU_{2}\textrm{~{}isospin}\big{>}$ states
$\big{|}\mathbf{8},\mathbf{1}\big{>}$ and
$\big{|}\mathbf{8},\mathbf{3}\big{>}$, formed out of
$\Lambda_{c}^{+}D_{s}^{-}$, $\Xi_{c}^{+}D^{-}$, $\Xi_{c}^{0}\bar{D}^{0}$ using
$SU_{3}$ isoscalar factors Kaeding:1995vq , are identical to those of
$\Lambda_{c}\bar{D}$, which is the $\big{|}\mathbf{8},\mathbf{2}\big{>}$ state
with the same spin structure. In a similar way, matrix elements involving
$\Lambda J/\psi$ and $\Sigma J/\psi$ ($\Lambda\eta_{c}$ and $\Sigma\eta_{c}$)
are the same as the corresponding terms in our previous paper involving
$NJ/\psi$ ($N\eta_{c}$). In summary, we have

$\displaystyle\big{<}\mathbf{8},\mathbf{1}\big{|}V\big{|}\mathbf{8},\mathbf{1}\big{>}$
$\displaystyle=\big{<}\mathbf{8},\mathbf{3}\big{|}V\big{|}\mathbf{8},\mathbf{3}\big{>}=A,$
(7) $\displaystyle\big{<}\mathbf{8},\mathbf{1}\big{|}V\big{|}\Lambda
J/\psi\big{>}$
$\displaystyle=\big{<}\mathbf{8},\mathbf{3}\big{|}V\big{|}\Sigma
J/\psi\big{>}=\frac{\sqrt{3}}{2}D,$ (8)
$\displaystyle\big{<}\mathbf{8},\mathbf{1}\big{|}V\big{|}\Lambda\eta_{c}\big{>}$
$\displaystyle=\big{<}\mathbf{8},\mathbf{3}\big{|}V\big{|}\Sigma\eta_{c}\big{>}=\frac{1}{2}D,$
(9)

where $A$ and $D$ are contact terms which are (somewhat) constrained by our
previous analysis of $\Lambda_{b}\to J/\psi\,pK^{-}$ decays Burns:2022uiv
(see below).

The $\Lambda_{c}^{+}D_{s}^{-}$, $\Xi_{c}^{+}D^{-}$, $\Xi_{c}^{0}\bar{D}^{0}$
basis states also combine into a flavor singlet, which implies an additional
contact term which is independent of the other contact terms Yan:2021nio . We
call this $A^{\prime}$, by analogy with $A$ above:

$\displaystyle\big{<}\mathbf{1},\mathbf{1}\big{|}V\big{|}\mathbf{1},\mathbf{1}\big{>}=A^{\prime}.$
(10)

Although the potentials are assumed to respect $SU_{3}$ and $SU_{2}$
symmetries, the transition amplitudes will not, because of mass differences
among the constituents with different flavors. Hence we formulate the
potential in the particle basis, and show the result in Table 1, where we have
introduced $\Delta=(A^{\prime}-A)/3$. Note that, as in our previous work, we
are assuming that potentials coupling two “closed-charm” states (such as
$\Lambda J/\psi\to\Lambda J/\psi$) are zero.

|  $\Lambda_{c}^{+}D_{s}^{-}$ |  $\Xi_{c}^{+}D^{-}$ |  $\Xi_{c}^{0}\bar{D}^{0}$ |  $\Lambda J/\psi$ |  $\Lambda\eta_{c}$ |  $\Sigma J/\psi$ |  $\Sigma\eta_{c}$
---|---|---|---|---|---|---|---
$\Lambda_{c}^{+}D_{s}^{-}$ |  $A+\Delta$ |  $\Delta$ |  $-\Delta$ |  $\frac{D}{\sqrt{2}}$ |  $\frac{D}{\sqrt{6}}$ |  0 |  0
$\Xi_{c}^{+}\bar{D}^{-}$ |  |  $A+\Delta$ |  $-\Delta$ |  $-\frac{D}{2\sqrt{2}}$ |  $-\frac{D}{2\sqrt{6}}$ |  $\frac{\sqrt{3}D}{2\sqrt{2}}$ |  $\frac{D}{2\sqrt{2}}$
$\Xi_{c}^{0}\bar{D}^{0}$ |  |  |  $A+\Delta$ |  $\frac{D}{2\sqrt{2}}$ |  $\frac{D}{2\sqrt{6}}$ |  $\frac{\sqrt{3}D}{2\sqrt{2}}$ |  $\frac{D}{2\sqrt{2}}$
$\Lambda J/\psi$ |  |  |  |  0 |  0 |  0 |  0
$\Lambda\eta_{c}$ |  |  |  |  |  0 |  0 |  0
$\Sigma J/\psi$ |  |  |  |  |  |  0 |  0
$\Sigma\eta_{c}$ |  |  |  |  |  |  |  0
Table 1: Contact terms in the $1/2^{-}$ channel.

With our model for the potentials, analysis of $\Lambda_{b}\to J/\psi\,pK^{-}$
decays Burns:2022uiv  indicates that $D$ is constrained very roughly to be a
number of order 1 GeV-2, and we will adopt that value in this work. The
contact term $A$ is not well-constrained by $\Lambda_{b}\to J/\psi\,pK^{-}$
decays, although we know that it cannot be large and negative, as it would
imply $\Lambda_{c}\bar{D}^{(*)}$ bound states, which are apparently not seen
in the data. Hence in this work we consider $A\geq 0$.

##  III Results

We first attempt a fit with only diagram $T_{1}$ (fixing $g_{2}=0$),
corresponding to the conventional molecular scenario. To understand the
possibilities, we note that in the isospin basis, the diagonal
$\Xi_{c}\bar{D}$ potentials are

$\displaystyle V(\Xi_{c}\bar{D},I=0)$ $\displaystyle=A+2\Delta,$ (11)
$\displaystyle V(\Xi_{c}\bar{D},I=1)$ $\displaystyle=A.$ (12)

Comparing to the $\Lambda_{c}^{+}D_{s}^{-}$ potential in Table 1, it suggests
there may be some region of parameter space (with $\Delta<0$) in which
$\Xi_{c}\bar{D}$ ($I=0$) binds, but not $\Lambda_{c}^{+}D_{s}^{-}$. But in
practice we find the opposite, namely a prominent signal at (or below)
$\Lambda_{c}^{+}D_{s}^{-}$ threshold, rather than $\Xi_{c}\bar{D}$. One reason
is that $\Delta$ not only contributes to the diagonal potential, but also the
off-diagonal coupling between the $\Lambda_{c}^{+}D_{s}^{-}$ and
$\Xi_{c}\bar{D}$ channels, generating an effective attraction in the lower
channel ($\Lambda_{c}^{+}D_{s}^{-}$), and also broadening any peak associated
with $\Xi_{c}\bar{D}$ (increasing its decay width). Ref Yan:2022wuz  also
observed that binding in $\Xi_{c}\bar{D}$ inevitably also implies binding in
$\Lambda_{c}^{+}D_{s}^{-}$. We also remark that, in our model, the
$\Lambda_{c}^{+}D_{s}^{-}$ state is always more prominent in the amplitude,
since it is produced directly in the triangle diagrams – in contrast to
$\Xi_{c}\bar{D}$, which arises in diagram $T_{1}$ only through final-state
interactions (via $\Delta$). All of this is problematic for the molecular
scenario.

Hence we proceed with the full model, including both diagrams $T_{1}$ and
$T_{2}$. The results are much better: in Fig. 2 we give two illustrative
examples of successful fits where, in both cases, the bulk of the distribution
is captured by diagram $T_{1}$, while the $P_{\psi s}^{\Lambda}(4338)$ peak is
from diagram $T_{2}$, due to the anticipated triangle singularity. The shape
of the triangle peak is sensitive to the chosen width of the virtual
$\bar{\Sigma}_{c}$, which in turn has implications for the contact terms $A$
and $\Delta$. The two parameter sets illustrated in Fig. 2 are summarized in
Table 2.

|  $\Gamma(\bar{\Sigma}_{c})$ / MeV |  $A$ / GeV-2 |  $\Delta$ / GeV-2
---|---|---|---
Set A |  70 |  6 | $-7$
Set B |  15 |  0 | $-1$
Table 2: Parameter sets A (green) and B (blue) in Fig 2.

The effect of a broader $\bar{\Sigma}_{c}$ is to broaden the triangle peak in
diagram $T_{2}$. Hence in Set A, in order to obtain a sufficiently sharp peak,
we adjust $A$ and $\Delta$ to effectively add some attraction to the
$\Xi_{c}\bar{D}$ channel, but in such a way that we do not have a bound state
in $\Xi_{c}\bar{D}$ or $\Lambda_{c}^{+}D_{s}^{-}$. With the narrower
$\bar{\Sigma}_{c}$ of Set B, the triangle peak is already sharp enough that no
extra attraction is really needed, although with our chosen $\Delta=-1$ GeV-2
there is a small attraction, as well as coupling between the $\Xi_{c}\bar{D}$
and $\Lambda_{c}^{+}D_{s}^{-}$ channels. Note that the values of $A$ and
$\Delta$ are not very well constrained by the fits.

To better understand our results, let us consider one of the fits (Set A) in
some more detail. The fit parameters are

$\displaystyle b$ $\displaystyle=0.93+6.63i\ \textrm{GeV}^{-2},$ (13)
$\displaystyle g_{1}$ $\displaystyle=5766,$ (14) $\displaystyle g_{2}$
$\displaystyle=-316.1.$ (15)

In Fig. 2 we show separately (with thin dotted lines) the contributions from
diagrams $T_{1}$ and $T_{2}$, and the background. Notice that $T_{1}$ nicely
captures the overall shape of the distribution, including the peak around the
$\Lambda_{c}^{+}D_{s}^{-}$ threshold which, as mentioned previously, is a
natural consequence of the assumed dominance of the color-favored triangle
diagram. While it is amusing that the peak corresponds to the single data
point near 140 candidates, we regard this as a fluke. It is also reassuring
that the background is relatively small compared to $T_{1}$, and this is as
expected, since tree-level production of the final state is color-suppressed.

We are also satisfied that the fit has $|g_{2}|<<|g_{1}|$, since it is
consistent with expectations that the electroweak vertex in $T_{2}$ is
suppressed compared to $T_{1}$ by at least the number of colors, with further
suppression expected due to the orbital excitation in $\bar{\Sigma}_{c}$. As
has been anticipated, despite a small coupling $g_{2}$, diagram $T_{2}$ still
makes a prominent contribution to the fit because of the enhancement due to
the triangle singularity.

The chi-squared value for the fit with Set A parameters is 1.9. This can be
improved substantially by using a more complicated background model; however,
our purpose is not to obtain a very precise fit, but to establish the physical
mechanism that explains the data.

Figure 2: The $J/\psi\Lambda$ invariant mass spectrum for parameter sets A
(green) and B (blue), compared to the experimental data from ref.
ChenParticleZoo . The thin dotted lines show the separate contributions of
$T_{1}$, background, and $T_{2}$ (top to bottom), for parameter set A.

From the matrix elements in Table 1, we can make some general observations
about the prospects of observing $P_{\psi s}^{\Lambda}(4338)$ in other final
states. (Detailed predictions for distributions in various final states are
difficult, without a model for the backgrounds in each case.)

For example, in $B^{-}\to\eta_{c}\,\Lambda\,\bar{p}$ decays, we expect a
$P_{\psi s}^{\Lambda}(4338)$ signal in the $\eta_{c}\,\Lambda$ distribution,
but suppressed in comparison to the $J/\psi\,\Lambda$ signal by a factor of
approximately 3 (in rate).

As noted previously, we expect $P_{\psi s}^{\Lambda}(4338)$ to exhibit isospin
mixing, as the $\Xi_{c}\bar{D}$ pair in diagram $T_{2}$ is a linear
combination of isospin 0 and 1. If we ignore the (small) additional
contribution to mixing due to the charged/neutral mass difference, and assume
that $P_{\psi s}^{\Lambda}(4338)$ is solely due to the triangle singularity
(in the sense that there are no non-perturbative final state interactions in
$\Xi_{c}\bar{D}$), we find that in diagram $T_{2}$, the signals in isospin 1
modes ($J/\psi\,\Sigma^{0}$ and $\eta_{c}\Sigma^{0}$) are suppressed compared
to the corresponding isospin 0 modes ($J/\psi\,\Lambda$ and
$\eta_{c}\Lambda$), but only by a factor of 3 in rate (before accounting for
phase space differences). Isospin mixing in this system is therefore a very
significant effect. Moreover, the $P_{\psi s}^{\Lambda}(4338)$ peak will
feature more prominently in the $J/\psi\,\Sigma^{0}$ and $\eta_{c}\Sigma^{0}$
distributions since (at leading order) they receive no contribution from the
diagram $T_{1}$, which dominates the $J/\psi\Lambda$ spectrum.

But the factor of three only applies in the perturbative limit, and indeed we
find that non-perturbative final state interactions can change the outcome
considerably. For example, in parameter set A (Table 2), the values of $A$ and
$\Delta$ imply attraction in isospin 0, but repulsion in isospin 1 – see
equations (11) and (12). So by introducing final state interactions to enhance
the triangle peak in isospin 0, we effectively suppress the peak in isospin 1.
Indeed we have verified that with our parameter set A, the suppression of the
isospin 1 peak is very much stronger than the factor of 3 which applies in the
perturbative limit.

Because the magnitude of isospin mixing is correlated with the parameters $A$
and $\Delta$, future experimental measurement of isospin 1 modes could be used
to constrain these parameters, which cannot be determined with current data.

##  IV Conclusion

We have demonstrated that a simple model based on the assumed dominance of
color-favored weak transitions, and approximate heavy-quark and $SU_{3}$
flavor symmetries, can describe the $B^{-}\to J/\psi\,\Lambda\,\bar{p}$ data,
including the prominent $P_{\psi s}^{\Lambda}(4338)$ peak. Notably, we are not
assuming a molecular nature for $P_{\psi s}^{\Lambda}(4338)$ – on the
contrary, we showed that the molecular scenario does not work. Our conclusions
(and methods) are very similar to those of our recent analysis of
$\Lambda_{b}\to J/\psi\,pK^{-}$ data Burns:2021jlu ; Burns:2022uiv , in which
we argued that the $\Sigma_{c}\bar{D}^{*}$ molecular interpretation of
$P_{c}(4457)$ is problematic, but that several viable alternatives give an
excellent fit to data. Our results underline the importance of exploring
alternatives to the prevailing molecular interpretation of states which are
located at thresholds.

###### Acknowledgements.

Swanson’s research was supported by the U.S. Department of Energy under
contract DE-SC0019232.

## References

  * (1) C. Chen and E. Spadaro Norella,  LHCb-PAPER-2022-031
https://indico.cern.ch/event/1176505/.

  * (2) M. Karliner and J. R. Rosner,  (2022), arXiv:2207.07581.
  * (3) F.-L. Wang and X. Liu,  (2022), arXiv:2207.10493.
  * (4) M.-J. Yan, F.-Z. Peng, M. S. Sánchez, and M. Pavon Valderrama,  (2022), arXiv:2207.11144.
  * (5) B. Wang, L. Meng, and S.-L. Zhu,  Phys. Rev. D 101, 034018 (2020), arXiv:1912.12592.
  * (6) C. W. Xiao, J. Nieves, and E. Oset,  Phys. Lett. B 799, 135051 (2019), arXiv:1906.09010.
  * (7) F.-Z. Peng, M.-Z. Liu, Y.-W. Pan, M. Sanchez Sanchez, and M. P. Valderrama,  (2019), arXiv:1907.05322.
  * (8) F.-Z. Peng, M.-J. Yan, M. Sánchez Sánchez, and M. P. Valderrama,  Eur. Phys. J. C 81, 666 (2021), arXiv:2011.01915.
  * (9) H.-X. Chen, W. Chen, X. Liu, and X.-H. Liu,  Eur. Phys. J. C 81, 409 (2021), arXiv:2011.01079.
  * (10) X.-K. Dong, F.-K. Guo, and B.-S. Zou,  Progr. Phys. 41, 65 (2021), arXiv:2101.01021.
  * (11) M.-J. Yan, F.-Z. Peng, M. S. Sánchez, and M. P. Valderrama,  Eur. Phys. J. C 82, 574 (2022), arXiv:2108.05306.
  * (12) K. Chen, R. Chen, L. Meng, B. Wang, and S.-L. Zhu,  Phys. Rev. D 103, 074007 (2021), arXiv:2109.13057.
  * (13) R. Chen,  Phys. Rev. D 103, 054007 (2021), arXiv:2011.07214.
  * (14) M.-Z. Liu, Y.-W. Pan, and L.-S. Geng,  Phys. Rev. D 103, 034003 (2021), arXiv:2011.07935.
  * (15) J.-T. Zhu, L.-Q. Song, and J. He,  Phys. Rev. D 103, 074007 (2021), arXiv:2101.12441.
  * (16) C. W. Xiao, J. J. Wu, and B. S. Zou,  Phys. Rev. D 103, 054016 (2021), arXiv:2102.02607.
  * (17) R. Chen and X. Liu,  Phys. Rev. D 105, 014029 (2022), arXiv:2201.07603.
  * (18) LHCb, R. Aaij et al.,  Sci. Bull. 66, 1278 (2021), arXiv:2012.10380.
  * (19) LHCb, R. Aaij et al.,  Phys. Rev. Lett. 115, 072001 (2015), arXiv:1507.03414.
  * (20) LHCb, R. Aaij et al.,  Phys. Rev. Lett. 122, 222001 (2019), arXiv:1904.03947.
  * (21) LHCb, R. Aaij et al.,  Phys. Rev. Lett. 117, 082002 (2016), arXiv:1604.05708.
  * (22) M.-Z. Liu et al.,  Phys. Rev. Lett. 122, 242001 (2019), arXiv:1903.11560.
  * (23) M.-Z. Liu et al.,  Phys. Rev. D 103, 054004 (2021), arXiv:1907.06093.
  * (24) M.-L. Du et al.,  Phys. Rev. Lett. 124, 072001 (2020), arXiv:1910.11846.
  * (25) S. Sakai, H.-J. Jing, and F.-K. Guo,  Phys. Rev. D 100, 074007 (2019), arXiv:1907.03414.
  * (26) M. P. Valderrama,  Phys. Rev. D 100, 094028 (2019), arXiv:1907.05294.
  * (27) C. W. Xiao, J. Nieves, and E. Oset,  Phys. Rev. D100, 014021 (2019), arXiv:1904.01296.
  * (28) F.-Z. Peng, M.-Z. Liu, M. Sánchez Sánchez, and M. Pavon Valderrama,  Phys. Rev. D 102, 114020 (2020), arXiv:2004.05658.
  * (29) M.-L. Du et al.,  JHEP 08, 157 (2021), arXiv:2102.07159.
  * (30) T. J. Burns and E. S. Swanson,  Eur. Phys. J. A 58, 68 (2022), arXiv:2112.11527.
  * (31) T. J. Burns and E. S. Swanson,  Phys. Rev. D100, 114033 (2019), arXiv:1908.03528.
  * (32) T. J. Burns and E. S. Swanson,  (2022), arXiv:2207.00511.
  * (33) T. J. Burns and E. S. Swanson,  Phys. Lett. B 813, 136057 (2021), arXiv:2008.12838.
  * (34) Belle, Y. B. Li et al.,  Phys. Rev. Lett. 122, 082001 (2019), arXiv:1811.09738.
  * (35) F.-K. Guo, X.-H. Liu, and S. Sakai,  Prog. Part. Nucl. Phys. 112, 103757 (2020), arXiv:1912.07030.
  * (36) Particle Data Group, P. A. Zyla et al.,  PTEP 2020, 083C01 (2020).
  * (37) L. Meng, B. Wang, and S.-L. Zhu,  (2022), arXiv:2208.03883.
  * (38) E. S. Swanson,  Phys. Lett. B588, 189 (2004), arXiv:hep-ph/0311229.
  * (39) F. E. Close and P. R. Page,  Phys. Lett. B578, 119 (2004), arXiv:hep-ph/0309253.
  * (40) T. J. Burns,  Eur. Phys. J. A51, 152 (2015), arXiv:1509.02460.
  * (41) F.-K. Guo, H.-J. Jing, U.-G. Meissner, and S. Sakai,  Phys. Rev. D99, 091501 (2019), arXiv:1903.11503.
  * (42) Z.-H. Guo and J. A. Oller,  Phys. Lett. B793, 144 (2019), arXiv:1904.00851.
  * (43) T. Burns and E. Swanson,  Phys. Rev. D 103, 014004 (2021), arXiv:2009.05352.
  * (44) L. A. Copley, N. Isgur, and G. Karl,  Phys. Rev. D20, 768 (1979),  [Erratum: Phys. Rev.D23,817(1981)].
  * (45) C. Chen, X.-L. Chen, X. Liu, W.-Z. Deng, and S.-L. Zhu,  Phys.Rev. D75, 094017 (2007), arXiv:0704.0075.
  * (46) H.-Y. Cheng and C.-K. Chua,  Phys. Rev. D75, 014006 (2007), arXiv:hep-ph/0610283.
  * (47) H.-Y. Cheng and C.-K. Chua,  Phys. Rev. D92, 074014 (2015), arXiv:1508.05653.
  * (48) K.-L. Wang and X.-H. Zhong,  Chin. Phys. C 46, 2 (2022), arXiv:2110.12443.
  * (49) M.-Z. Liu, F.-Z. Peng, M. Sanchez Sanchez, and M. P. Valderrama,  Phys. Rev. D98, 114030 (2018), arXiv:1811.03992.
  * (50) T. A. Kaeding,  (1995), arXiv:nucl-th/9502037.


2020q2 |  | 1.9359 |  | 1.6977 |  | 1.9417 |  | 1.8175 |  | 1.9417
2020q3 |  | 1.6059 |  | 1.5965 |  | 1.4484 |  | 1.7514 |  | 1.4484
2020q4 |  | 0.9909 |  | 1.6540 |  | 1.9467 |  | 1.7858 |  | 1.9467
2021q1 |  | –0.1663 |  | 1.7455 |  | 1.7464 |  | 1.8576 |  | 1.7464
2021q2 |  | –4.6085 |  | 1.7584 |  | 0.4562 |  | 2.0774 |  | 0.4562
2021q3 |  | 8.6921 |  | 2.2141 |  | 1.3667 |  | 2.5984 |  | 1.3667
2021q4 |  | 1.5211 |  | 2.0870 |  | 0.7506 |  | 2.4724 |  | 0.7506
2022q1 |  | 1.8425 |  | 2.2670 |  | 1.5652 |  | 2.5388 |  | 1.5652
2022q2 |  | 2.8290 |  | 8.3146 |  | 2.4832 |  | 3.6200 |  | 2.4832
2022q3 |  | 0.6644 |  | 3.5159 |  | 2.7159 |  | 2.9876 |  | 2.7159
2022q4 |  | 2.1661 |  | 4.2029 |  | 3.7504 |  | 3.1511 |  | 3.7504
2023q1 |  | 3.1894 |  | 4.8976 |  | 4.6815 |  | 4.0785 |  | 4.7542
2023q2 |  | 4.9977 |  | 4.5467 |  | 5.7679 |  | 4.4758 |  | 5.7679
2023q3 |  | 5.4954 |  | 4.6832 |  | 5.6688 |  | 4.6143 |  | 5.6688
2023q4 |  | 4.9713 |  | 4.7016 |  | 5.2734 |  | 4.9634 |  | 5.2734
2024q1 |  | 2.9576 |  | 3.6366 |  | 3.1912 |  | 4.2441 |  | 3.1912
RMSFE |  | 4.5021 |  | 2.4643 |  | 3.2470 |  | 2.9795 |  | 3.2402

___Notes_ : The ”Actual” column gives the realized inflation rate at period
$t,$ the other columns the forecasts for each model. RMSFE is the Root Mean
Square Forecast Error. Lasso models are estimated using the code of Chudik,
Kapetanios, and Pesaran (2018).

Figure S.3: Four Quarters Ahead Forecast Plots

## References

Chudik, Kapetanios, and Pesaran (2018) “A One Covariate at a Time, Multiple
Testing Approach to Variable Selection in High-Dimensional Linear Regression
Models,”  _Econometrica_ , Vol. 86, No. 4, 1479-1512.


Induced Gravitational Waves from Non-attractor Inflation
and NANOGrav data

Amin Nassiri-Rad<EMAIL_ADDRESS>Kosar Asadi<EMAIL_ADDRESS>
School of Astronomy, Institute for Research in Fundamental Sciences (IPM)
P. O. Box 19395-5531, Tehran, Iran

###### Abstract

In this paper, we investigate the scalar-induced gravitational waves in
single-field non-attractor inflation for the Pulsar Timing Arrays data. Our
model comprises three phases of inflation: the first and third phases are
slow-roll inflation, while the second phase is a period of non-attractor
inflation. We analyze the model’s predictions for various values of the sound
speed $c_{s}$ and examine the sharp transitions to the final attractor phase.
Furthermore, we study the model’s predictions for NANOGrav observations and
future gravitational wave observations. We also calculate the non-Gaussianity
parameter $f_{NL}$ for the non-attractor setup with a general sound speed and
the sharpness parameter.

##  1 Introduction

The detection of gravitational waves (GWs) has indeed revolutionized our
understanding of physics on large scales and in non-linear regimes. Recent
advancements have led to significant progress in detecting stochastic GWs in
the range of 10 nHz, as observed by NANOGrav [1] and other experiments [2, 3,
4]. The NANOGrav data provides valuable insights into the astrophysical
origins of GWs, including the merging of binary black holes. Another
intriguing possibility is the existence of induced gravitational waves, which
are sourced by scalar perturbations on small scales during inflation [14, 13,
12, 11, 10, 9, 8, 7, 6, 5].

To explain the NANOGrav data via cosmological sources, one needs to amplify
the scalar perturbations during inflation [16, 17, 18, 19, 20, 21, 22, 23, 24,
25]. The amplification of scalar perturbations is also a mechanism for the
production of primordial black holes (PBHs), which may explain dark matter as
well[26, 27, 28, 29, 30]. The GW perturbations are then sourced by the square
of the scalar perturbations. If the amplitude of the scalar perturbations at
small scales is much higher than the corresponding amplitude at CMB scales,
then there is a chance that the induced GWs are large enough to explain
NANOGrav data.

The Ultra slow-roll (USR) inflation model is a well-known approach to amplify
scalar perturbations, which can lead to the formation of PBHs [31, 32, 33, 34,
35, 36, 37, 38, 39]. In this model, the potential is constant, resulting in an
exponential decrease in the velocity of the inflaton over time. Consequently,
this leads to an increase in the curvature perturbation [40, 41, 42, 43]. It
is important to note that the USR model violates the Maldacena consistency
condition [44, 45]. Several studies have explored this violation, including
[47, 48, 49, 51, 52, 53, 54, 55]. According to [43], the non-Gaussianity
parameter in this model is $f_{NL}=\frac{5}{2}$. In [46] this result is
extended by introducing a sharpness parameter that characterizes the final
attractor phase.

In an alternative scenario, the perturbations can experience growth during the
non-attractor phase, even with an arbitrary sound speed $c_{s}$ [48, 49]. In
this study, we extend the approach proposed in [56] by considering the
arbitrary sound speed $c_{s}$ as an additional degree of freedom. Our
objective is to investigate how this model can account for the NANOGrav data.
In our model, the Lagrangian of the inflaton is an arbitrary function of the
kinetic term $X$, denoted as $P(X)$, with a constant potential. The sound
speed can have non-trivial effects on the generalization of the results
obtained in [56]. Specifically, our model comprises three inflationary phases:
two slow-roll phases at the beginning and end, and a non-attractor phase in
between. During the slow-roll phases, the sound speed is set equal to unity.
In this work, we also explore non-Gaussianity and demonstrate that it depends
on both the sharpness parameter $h$ and the number of e-folds associated with
the non-attractor phase. Notably, when $c_{s}=1$, our results reduce to the
standard result obtained in [46]. Furthermore, if the transition to the final
attractor phase is very mild, non-Gaussianity approaches zero.

The researches on ultra slow-roll inflation has recently focused on the loop
effects of this model on the CMB scale [57, 58, 59, 60, 61, 62, 63, 64, 65,
66, 67, 68, 69, 70]. The main message conveyed by these works is that loop
corrections can become significant and affect long CMB scales perturbations,
depending on the duration of ultra slow-roll inflation. In particular, it is
claimed in [63] that if the transition to the final attractor phase is mild,
then loop corrections will become harmless. However, other criticisms have
been raised that loop corrections are absent [59, 60, 69]. This question
remains open, and future studies will be needed to investigate the loop
effects in our model.

The paper is structured as follows. In Section 2, we present our model and
generalize the results of [56] to the case with a general sound speed. We
calculate the non-Gaussianity for this model and study its consistency in
confrontation with the observations of PBHs. In Section 3, we compare our
results with the NANOGrav data and discuss the differences compared to the USR
model. Our main focus is to investigate the predictions of the model for
various values of $c_{s}$ and see which ones fit better to the data
qualitatively. We conclude in Section 4.

##  2 The Model

In this section, we introduce our model for the scalar-induced gravitational
waves as the source of pulsar timing arrays (PTAs) observations. Our setup
comprises three stages: two slow-roll inflation stages at the beginning and
end, and a non-attractor phase in the middle. This model is a generalization
of the work presented in [56], where the middle phase allows for an arbitrary
sound speed, $c_{s}$. Specifically, the Lagrangian in the middle phase is
given by

$\mathcal{L}_{2}=P(X)-V_{0},$ (2.1)

where $V_{0}$ is a constant, and $X=\frac{\dot{\phi}^{2}}{2}$ is the kinetic
energy. The sound speed in terms of $P(X)$ is given by [50]

$c_{s}^{2}\equiv\frac{P_{,X}}{P_{,X}+2XP_{,XX}}.$ (2.2)

As $P(X)=X$ in the first and third stages, the sound speed is set to unity
($c_{s}=1$) in these phases. In the non-attractor phase, the evolution of
$\phi$ at the background level is given by [49]

$\ddot{\phi}+3H\dot{\phi}c_{s}^{2}=0.$ (2.3)

In the first slow-roll phase, the mode function in terms of conformal time
$\tau$ is given by

$\mathcal{R}_{1}=\frac{He^{-ik\tau}(1+ik\tau)}{2\sqrt{k^{3}\text{$\epsilon_{i}$}}}.$
(2.4)

Here we assume that $\epsilon_{i}$ represents the slow roll parameter at the
first phase, which is assumed to be nearly constant. We also consider the
Bunch-Davies initial condition for the primary mode function. Throughout this
paper, we assume $M_{P}=1$ without loss of generality. In the middle stage, as
the sound speed is arbitrary, the mode function will be more complicated.
Following [75], the mode function in the non-attractor phase evolves as

$v_{k}^{\prime\prime}+(c_{s}^{2}k^{2}-\frac{(\eta-2s+3)^{2}-1}{4\tau^{2}})v_{k}=0,$
(2.5)

where we define the variable $v_{k}=z\mathcal{R}_{k}$ , $z=\frac{2\epsilon
a^{2}}{c_{s}^{2}}$ and $s=\frac{\dot{c_{s}}}{Hc_{s}}$. Here, $a$ represents
the scale factor. Additionally, we have the first and second slow-roll
parameters, $\epsilon$ and $\eta$, respectively. They are defined as
$\epsilon=-\frac{\dot{H}}{H^{2}}$ and $\eta=\frac{\dot{\epsilon}}{H\epsilon}$
(a ’dot’ represents the differentiation with respect to the cosmic time). In
the non-attractor phase, we assume that the phase begins at $\tau=\tau_{i}$
and ends at $\tau=\tau_{e}$. Based on this assumption, it can be shown that
the first and second slow-roll parameters during this phase with a constant
potential are expressed as

$\epsilon_{n}(\tau)=\text{$\epsilon_{i}$}\left(\frac{\tau}{\tau_{i}}\right)^{-\eta},$
(2.6) $\eta=-3(1+c_{s}^{2}).$ (2.7)

The general solution to (2.5) in the non-attractor phase is given as

$\mathcal{R}_{2}=\frac{iHc_{s}\sqrt{-\tau^{3}}}{2\sqrt{\epsilon}_{i}}\left(\frac{\tau_{i}}{\tau}\right)^{-\frac{\eta}{2}}\bigg{[}\alpha_{k}H_{v}^{(1)}(-c_{s}k\tau)+\beta_{k}H_{v}^{(2)}(-c_{s}k\tau)\bigg{]},$
(2.8)

where $\nu=\frac{-3c_{s}^{2}}{2}$. $\alpha_{k}$ and $\beta_{k}$ are constants
that should be determined by the matching condition at $\tau=\tau_{i}$. As the
sound speed is different in the middle stage the matching conditions are as
follows[76]

$\begin{split}&R_{k}(\tau_{i}^{-})=R_{k}(\tau_{i}^{+})\\\
&R^{\prime}_{k}(\tau_{i}^{-})=\frac{R^{\prime}_{k}(\tau_{i}^{+})}{c_{s}^{2}}.\end{split}$
(2.9)

we can determine the values of $\alpha_{k}$ and $\beta_{k}$ as

$\alpha_{k}=\frac{ie^{-ik\tau_{i}}}{2\sqrt{-k^{3}\tau_{i}^{3}}}\bigg{[}k\tau_{i}\left(k\tau_{i}-i\right)H_{1+\nu}^{(2)}\left(-kc_{s}\tau_{i}\right)+ic_{s}\left(3+k\tau_{i}\left(k\tau_{i}+3i\right)\right)H_{\nu}^{(2)}\left(-kc_{s}\tau_{i}\right)\bigg{]},$
(2.10)
$\beta_{k}=-\frac{ie^{-ik\tau_{i}}}{2\sqrt{-k^{3}\tau_{i}^{3}}}\bigg{[}k\tau_{i}\left(k\tau_{i}-i\right)H_{1+\nu}^{(1)}\left(-c_{s}k\tau_{i}\right)+ic_{s}\left(3+k\tau_{i}\left(k\tau_{i}+3i\right)\right)H_{\nu}^{(1)}\left(-c_{s}k\tau_{i}\right)\bigg{]}.$
(2.11)

One can determine the mode function at the third stage by utilizing
$\alpha_{k}$ and $\beta_{k}$. It is possible to verify that when $c_{s}=1$,
the aforementioned expressions reduce to the findings of [63]. The first slow-
roll parameter in the final stage is given by [46]

$\epsilon(\tau)=\epsilon_{e}\left[\frac{h}{6}-\left(\frac{h}{6}+1\right)\frac{\tau^{3}}{\tau_{e}^{3}}\right]^{2}.$
(2.12)

Here, the parameter $h$ determines the sharpness of the transition to the
final attractor phase. A higher value of $|h|$ leads to a quicker transition.
In our model, $h$ is defined as[46, 63]

$h=-6\sqrt{\frac{\epsilon_{V}}{\epsilon_{e}}},$ (2.13)

where $\epsilon_{e}$ represents the slow-roll parameter at the end of the non-
attractor phase, and
$\epsilon_{V}=\frac{1}{2}(\frac{V^{\prime}(\phi_{e})}{V(\phi_{e})})^{2}$.

One can use these results to construct the mode function. The general mode
function in the final attractor phase is expressed as

$\mathcal{R}_{3}=\frac{H}{2\sqrt{k^{3}\epsilon(\tau)}}\bigg{[}c_{k}e^{-ik\tau}(1+ik\tau)+d_{k}e^{ik\tau}(1-ik\tau)\bigg{]}$
(2.14)

The constants $c_{k}$ and $d_{k}$ are determined by matching conditions
similar to (2.9) at $\tau=\tau_{e}$. Due to their complexity, we do not
provide explicit expressions for them. The power spectrum at the end of
inflation can be obtained using the following relation

$\mathcal{P}_{\mathcal{R}}=\frac{k^{3}}{2\pi^{2}}|\mathcal{R}_{3}|^{2}\Big{|}_{\tau=0}$
(2.15)

The two graphs in Fig. 1 illustrate the power spectrum for various
combinations of $c_{s}$ and $h$. To calculate the power spectrum, we select
$\tau_{i}$ such that the non-attractor phase commences about 16 e-folds after
the CMB mode exits the horizon. Additionally, we employ COBE normalization to
establish the power spectrum on CMB scales. In simpler terms, by satisfying
the matching condition, it can be demonstrated that on super-horizon scales,
we observe

$\lim_{k\rightarrow
0}\frac{k^{3}}{2\pi^{2}}|\mathcal{R}_{3}|^{2}=\frac{H^{2}}{8\pi^{2}\epsilon_{i}}$
(2.16)

The values of $H$ and $\epsilon_{i}$ are carefully selected such that equation
(2.16) yields a power spectrum $\mathcal{P}_{\mathcal{R}}\simeq 2\times
10^{-9}$ for $k_{CMB}=0.05\text{Mpc}^{-1}$. As depicted in Fig. 1, the power
spectrum exhibits similar characteristics to those described in [56], which
will be elaborated upon in detail. Let $\Delta N$ denote the e-fold duration
of the non-attractor phase. For $|\tau_{e}|\ll|\tau_{i}|$ and $k|\tau_{e}|\ll
1$ , it can be demonstrated that

$\begin{split}\mathcal{P}_{\mathcal{R}}=&\mathcal{P}_{\mathcal{\text{CMB}}}\frac{(6-h)^{2}}{h^{2}}e^{6c_{s}^{2}\Delta
N}\times\\\
&\Bigg{|}(1+ix)\,_{0}F_{1}\left(;\frac{3c_{s}^{2}}{2};-\frac{1}{4}x^{2}c_{s}^{2}\right)-\frac{1}{3}(3+x(x+3i))\,_{0}F_{1}\left(;\frac{3c_{s}^{2}}{2}+1;-\frac{1}{4}x^{2}c_{s}^{2}\right)\Bigg{|}^{2},\end{split}$
(2.17)

where $x=-k\tau_{i}$ and ${}_{0}F_{1}$ is the hypergeometric function. From
the above equation, we observe that the exponential growth of the power
spectrum is modified to $\exp(6c_{s}^{2}\Delta N)$. This indicates that the
peak of the non-attractor phase is shorter compared to the case where
$c_{s}=1$. Additionally, the power spectrum is proportional to
$\frac{(6-h)^{2}}{h^{2}}$. It is worth noting that the scaling behavior of the
power spectrum in the attractor phase is identical to that observed in ultra
slow-roll inflation. Another interesting feature to study is the behavior of
the power spectrum in the limit $|k\tau_{i}|\ll 1$, which pertains to modes
that exited the horizon at early times. Given the power spectrum at the end of
inflation, one can observe that the power spectrum at CMB scales is modified
by

$\mathcal{P}_{\mathcal{R}}(k,\tau=0)\simeq\mathcal{P}_{\text{CMB}}\bigg{[}1+\frac{2x^{2}\left(c_{s}^{2}+1\right)\left(6-h\right)e^{3c_{s}^{2}\Delta
N}}{h\left(3c_{s}^{2}+2\right)}\bigg{]}.$ (2.18)

As we see in Fig. 1 in the dip, the power spectrum decreases considerably. To
estimate the position of the dip, we expect Eq. (2.18) goes to zero at dip and
one has (2.18)

$k|\tau_{i}|\sim e^{-3c_{s}^{2}\Delta
N/2}\sqrt{\frac{h\left(3c_{s}^{2}+2\right)}{2\left(c_{s}^{2}+1\right)\left(h-6\right)}}.$
(2.19)

The above expression also demonstrates that the position of the dip shifts to
the right as the absolute value of $|h|$ increases.

Figure 1:  Left: The power spectrum for different values of $c_{s}$ and
$h=-6$. Other parameters are held constant to ensure that the peak for
different sound speeds does not grow significantly. Right: The power spectrum
for $c_{s}=\frac{3}{4}$ and different values of $h$.

One can also check the behaviour of power spectrum for the modes very deep
inside horizon. For $1\ll k|\tau_{e}|$ one can check that

$\begin{split}\mathcal{P}_{\mathcal{R}}\simeq&\frac{9H^{2}e^{3\Delta
N(1+c_{s}^{2})}}{8\pi^{2}h^{2}c_{s}^{2}\epsilon_{i}}\times\\\
&\Bigg{|}\left(c_{s}-1\right)e^{2ikc_{s}(\tau_{i}-\tau_{e})}\left(\cos\left(k\tau_{e}\right)+ic_{s}\sin\left(k\tau_{e}\right)\right)+\left(c_{s}+1\right)\left(\cos\left(k\tau_{e}\right)-ic_{s}\sin\left(k\tau_{e}\right)\right)\Bigg{|}^{2}\,.\end{split}$
(2.20)

The given expression oscillates for large values of k when $c_{s}\neq 1$.
However, if we consider $c_{s}=1$, the power spectrum remains constant for
$1\ll k|\tau_{e}|$, as depicted in the left panel of Figure 1. In particular,
the expression reduces to the following one by setting $c_{s}=1$

$\mathcal{P}_{\mathcal{R}}\simeq\frac{9H^{2}e^{6\Delta
N}}{2\pi^{2}h^{2}\epsilon_{i}}\,.$ (2.21)

To fit the parameters of the model with PTAs data, we also consider the
observed values for PBH formation, as shown in Fig. 2. To begin, we focus on
the mass fraction of PBHs, which is given by[71]

$\beta\simeq\frac{1}{2}\text{Erfc}\bigg{(}\frac{\mathcal{R}_{c}}{\sqrt{2\mathcal{P}_{\mathcal{R}}}}\bigg{)}.$
(2.22)

Here, $\mathcal{R}_{c}$ is an $O(1)$ parameter that denotes the critical value
of the curvature perturbation for PBH formation. In this model, we have chosen
$\mathcal{R}_{c}=1.59$. Additionally, the ratio of PBHs to the dark matter
density, $f_{PBH}$, is expressed as [28]

$f_{\text{PBH}}(M_{\text{PBH}})\simeq 2.7\times
10^{8}\bigg{(}\frac{M_{\text{PBH}}}{M_{\odot}}\bigg{)}^{-\frac{1}{2}}\beta.$
(2.23)

To determine $\frac{M_{\text{PBH}}}{M_{\odot}}$ in terms of the number of
e-folds, one can use the following relation [29]

$\frac{M_{\text{PBH}}}{M_{\odot}}=7.7\times
10^{17}\exp(-2(N_{\text{max}}-N_{\text{CMB}})).$ (2.24)

In our model, $N_{\text{max}}$ represents the number of e-folds at which the
power spectrum reaches its maximum and $N_{CMB}$ is the number of e-folds at
which CMB mode exits the horizon. For each value of $c_{s}$ and $h$, our model
has three free parameters: $h$, $\tau_{i}$ and $\Delta N$. We have carefully
selected these parameters such that for each $c_{s}$ the frequency of the peak
in the power spectrum is comparable to the frequency of PTAs data. In Fig. 2
we have presented constraints on $f_{\text{PBH}}$ for different values of
$c_{s}$. As the constraints on $f_{\text{PBH}}$ for different parameters were
somewhat similar, we have consolidated them into a single figure. When $h=-1$,
the amplitude of oscillations does not decay for large momentum values, as
observed in Figure 1. Since the approximation of $\beta$ in Equation (2.22)
may not be valid in this limit, we do not study the effect of mild
transitions. Before ending this section, we delve into the study of the
bispectrum for non-attractor models with arbitrary sound speed $c_{s}$ in the
next subsection.

Figure 2: This figure illustrates the behavior of the PBH mass fraction
($f_{\text{PBH}}$) for $c_{s}=\frac{3}{4}$ and $h=-12$. It is worth noting
that this behavior is similar for other values of $c_{s}$ and $h$, so we have
consolidated them into a single figure.

###  2.1 Bispectrum

In this subsection, we evaluate the bispectrum for the modes that exit the
horizon during the intermediate non-attractor phase. To calculate the
bispectrum, one can adopt a two-phase inflationary framework, where the first
phase corresponds to the non-attractor phase, and the second phase corresponds
to the slow-roll phase [63]. Since the mode functions during non-attractor
inflation do not exhibit scale invariance, we can consider the following
option for the mode function

$\mathcal{R}_{1}(\tau)=\frac{Hc_{s}\sqrt{-\pi\tau^{3}}}{2\sqrt{2\epsilon_{n}(\tau)}}\bigg{[}c_{+}H_{\nu}^{(1)}\left(-c_{s}k\tau
c_{s}\right)+c_{-}H_{\nu}^{(1)}\left(-c_{s}k\tau\right)\bigg{]}.$ (2.25)

The mode function in the second phase can be obtained using the following
expressions for the coefficients

$\begin{split}&c_{\pm}=\frac{1}{2}\left(\frac{k}{\lambda}\right)^{-\frac{3}{2}\left(c_{s}^{2}-1\right)}\left[\left(\frac{k}{\lambda}\right)^{3(c_{s}^{2}-1)}\pm
1\right]\,,\end{split}$ (2.26)

where $\lambda$ is an arbitrary conststant. It is straightforward to verify
that with these choices of coefficients, the mode function satisfies the
normalization condition

$c_{+}^{2}-c_{-}^{2}=1.$ (2.27)

The mode function in the second phase is then given by

$\mathcal{R}_{2}(\tau)=\frac{H}{2\sqrt{k^{3}\epsilon(\tau)}}\bigg{[}a_{k}e^{-ik\tau}(1+ik\tau)+b_{k}e^{ik\tau}(1-ik\tau)\bigg{]},$
(2.28)

with $a_{k}$ and $b_{k}$ to be determined by the matching conditions. Based on
the definition of $f_{NL}$, we have the following relation

$\left<\mathcal{R}_{k_{1}}\mathcal{R}_{k_{2}}\mathcal{R}_{k_{3}}\right>=(2\pi)^{3}\delta^{3}\bigg{(}\sum_{i}\overrightarrow{k_{i}}\bigg{)}\frac{6}{5}f_{NL}\big{(}P_{k_{1}}P_{k_{2}}+P_{k_{2}}P_{k_{3}}+P_{k_{3}}P_{k_{1}}\big{)},$
(2.29)

where $P_{k}=|R_{k}(\tau=0)|^{2}$. In order to calculate the bispectrum, we
require the Hamiltonian up to third order. Additionally, it is more convenient
to compute $f_{\text{NL}}$ using the effective field theory approach. It can
be shown that The third-order Hamiltonian, which neglects sub-leading gradient
interactions, can be expressed as follows [75]

$\begin{split}&H_{3}^{(1)}=-\frac{a^{2}\eta H^{3}\epsilon\pi^{\prime
2}\pi}{c_{s}(\tau)^{2}}\,,\\\
&H_{3}^{(2)}=\left[1-\frac{1}{c_{s}(\tau)^{2}}\right]aH^{2}\epsilon\pi^{\prime
3}\,,\\\
&H_{3}^{(3)}=2H^{2}a\epsilon\frac{c_{s}^{\prime}(\tau)}{c_{s}^{3}(\tau)}\pi^{\prime
2}\pi,\end{split}$ (2.30)

where $c_{s}(\tau)$ is the time dependent sound speed. The total third order
Hamiltonian is expressed as $H_{3}=H_{3}^{(1)}+H_{3}^{(2)}+H_{3}^{(3)}$. Here,
the prime denotes the derivative with respect to conformal time, and $\pi$
represents the Goldstone boson, which is related to the curvature perturbation
according to the relations established in [44, 73]

$\mathcal{R}_{k}=-H\pi_{k}+H\pi_{k}\dot{\pi}_{k}+\frac{\dot{H}}{2}\pi_{k}^{2}\,.$
(2.31)

Since our focus is on the three-point correlation function at the end of
inflation, the higher-order terms are suppressed within the slow-roll
approximation. Therefore, we can utilize the aforementioned relation at the
linear order, leading to an approximate expression for the three-point
correlation function as follows

$\left<\mathcal{R}_{k_{1}}\mathcal{R}_{k_{2}}\mathcal{R}_{k_{3}}\right>\simeq-H^{3}\left<\pi_{k_{1}}\pi_{k_{2}}\pi_{k_{3}}\right>.$
(2.32)

By employing the in-in formalism, we can express the three-point correlation
function as follows

$\left<\pi_{k_{1}}(0)\pi_{k_{2}}(0)\pi_{k_{3}}(0)\right>=-2\imaginary\int_{\tau_{i}}^{0}d\tau\left<H_{3}(\tau)\pi_{k_{1}}(0)\pi_{k_{2}}(0)\pi_{k_{3}}(0)\right>.$
(2.33)

The integral mentioned above can be divided into three parts:
$\tau_{i}<\tau<\tau_{e}$, $\tau=\tau_{e}$ and$\tau_{e}<\tau<0$. At
$\tau=\tau_{e}$ the Hamiltonian at the transition time contributes to the
three-point function. To compute the contribution of each part, we utilize the
mode functions corresponding to the non-attractor phase in the first part, and
those associated with the slow-roll phase in the second part. For
$\tau_{i}<\tau<\tau_{e}$ we have the non-attractor phase and
$c_{s}(\tau)=c_{s}$. In this part both $H_{3}^{(1)}$ and $H_{3}^{(2)}$
contribute to the three-point function and we have

$\frac{\left<\mathcal{R}_{k_{1}}\mathcal{R}_{k_{2}}\mathcal{R}_{k_{3}}\right>^{(1)}_{\text{non-
att}}}{2\pi^{2}\mathcal{P}_{\mathcal{R}}}=(2\pi)^{3}\delta^{3}\bigg{(}\sum_{i}\overrightarrow{k_{i}}\bigg{)}\frac{-3h\left(c_{s}^{2}+1\right)(h+12)}{2(h-6)^{2}}\left(\frac{1}{k_{1}^{3}}\frac{1}{k_{2}^{3}}+2\text{
perms}\right),$ (2.34)

$\frac{\left<\mathcal{R}_{k_{1}}\mathcal{R}_{k_{2}}\mathcal{R}_{k_{3}}\right>^{(2)}_{\text{non-
att}}}{2\pi^{2}\mathcal{P}_{\mathcal{R}}}=(2\pi)^{3}\delta^{3}\bigg{(}\sum_{i}\overrightarrow{k_{i}}\bigg{)}\frac{9h^{2}c_{s}^{2}\left(c_{s}^{2}-1\right)}{4(h-6)^{2}}\left(\frac{1}{k_{1}^{3}}\frac{1}{k_{2}^{3}}+2\text{
perms}\right),$ (2.35)

with $\mathcal{P}_{\mathcal{R}}=\lim_{k\rightarrow
0}\frac{k^{3}}{2\pi^{2}}|\mathcal{R}_{3}|^{2}$. To clarify, we use the term
‘non-att’ to refer to the non-attractor phase. On the other hand for
$\tau_{e}<\tau<0$ we have the slow roll phase. In this stage only
$H_{3}^{(1)}$ contributes to the three-point function with $c_{s}(\tau)=1$ .
Computing the three-point function in this part yields

$\frac{\left<\mathcal{R}_{k_{1}}\mathcal{R}_{k_{2}}\mathcal{R}_{k_{3}}\right>^{(1)}_{\text{SR}}}{2\pi^{2}\mathcal{P}_{\mathcal{R}}}=(2\pi)^{3}\delta^{3}\bigg{(}\sum_{i}\overrightarrow{k_{i}}\bigg{)}\frac{6h(h+6)e^{3\Delta
N\left(c_{s}^{2}-1\right)}}{(h-6)^{2}}\left(\frac{1}{k_{1}^{3}}\frac{1}{k_{2}^{3}}+2\text{
perms}\right).$ (2.36)

Finally, we calculate the contribution of the transition term. Since the sound
speed changes during the transition, we can express $c_{s}(\tau)^{2}$ as

$c_{s}(\tau)^{2}=c_{s}^{2}\theta(\tau_{e}-\tau)+\theta(\tau-\tau_{e}),$ (2.37)

with $\theta(x)$ being the step function. Writing
$H_{3}^{(3)}=H^{2}a\epsilon\frac{(c_{s}(\tau)^{2})^{\prime}}{c_{s}^{4}(\tau)}\pi^{\prime
2}\pi$ and using (2.9) at $\tau=\tau_{e}$, the above equation gives the
contribution of transition term to the three-point function as

$\frac{\left<\mathcal{R}_{k_{1}}\mathcal{R}_{k_{2}}\mathcal{R}_{k_{3}}\right>^{(3)}_{\text{transition}}}{2\pi^{2}\mathcal{P}_{\mathcal{R}}}=(2\pi)^{3}\delta^{3}\bigg{(}\sum_{i}\overrightarrow{k_{i}}\bigg{)}\frac{3h(2h+6)}{(h-6)^{2}}e^{-3\Delta
N(1-c_{s}^{2})}(c_{s}^{2}-1)\left(\frac{1}{k_{1}^{3}}\frac{1}{k_{2}^{3}}+2\text{
perms}\right).$ (2.38)

Taking into account all the above contributions the non-Gaussianity then reads
as

$f_{NL}=\frac{5h}{2(h-6)^{2}}\left[2\left((h+3)c_{s}^{2}+3\right)e^{-3\Delta
N\left(1-c_{s}^{2}\right)}+\frac{\left(-3hc_{s}^{4}+(h-24)c_{s}^{2}-2(h+12)\right)}{4}\right]\,.$
(2.39)

The expression mentioned above violates the Maldacena consistency condition,
which is similar to what happens in ultra slow-roll inflation. It is
noteworthy that when $c_{s}=1$, the expression reduces to the result obtained
in [63]. Furthermore, we can analyze the sharp transition limit where
$h\rightarrow-\infty$. In this limit, the expression simplifies to

$f_{NL}=\frac{5}{8}\left[c_{s}^{2}\left(8e^{-3\Delta
N\left(1-c_{s}^{2}\right)}+1\right)-3c_{s}^{4}-2\right]\,.$ (2.40)

It should be noted that the aforementioned outcome should not be mistaken for
the outcome obtained in [49]. In that study, the non-attractor phase is
assumed to have $\eta=-6$, which results in a non-constant potential. In
contrast, in the present analysis, we observe that for $|h|\ll 1$, $f_{NL}$
approaches zero, which is similar to the result obtained in the case of ultra
slow-roll inflation [63].

According to [72], the Gaussian bound for the PBH constraints is applicable
when $|f_{NL}|\mathcal{P}_{\mathcal{R}}\ll 1$. For instance, with $h=-6$,
$c_{s}=\frac{1}{4}$, and $\Delta N=3$, we obtain $f_{NL}\simeq 0.3$, and the
Gaussian approximation constraint is satisfied with
$\mathcal{P}_{\mathcal{R}}=\mathcal{O}(0.01)$. It is worth noting that the
non-Gaussianity parameter presented in (2.39) corresponds to a model with two
phases. However, if one were to inquire about the non-Gaussianity parameter in
a model with three phases, our analysis for this case confirms that
$f_{NL}=\mathcal{O}(\epsilon)$, which is negligible. Therefore, the
aforementioned criterion is satisfied once again.

##  3 PTAs observations and the scalar induced gravitational waves

The evolution of gravitational wave perturbations follows the same description
as presented in [74]. Specifically, the evolution of the Fourier mode of
gravitational waves can be described as follows

$h^{\lambda^{\prime\prime}}_{k}(\tau)+2\mathcal{H}h^{\lambda^{\prime}}_{k}(\tau)+k^{2}h^{\lambda}_{k}(\tau)=4S_{k}^{\lambda}(\tau),$
(3.1)

Here, $\mathcal{H}=\frac{a^{\prime}(\tau)}{a(\tau)}$ represents the conformal
Hubble parameter, and $S_{k}^{\lambda}(\tau)$ corresponds to the source term,
which is given by

$S_{k}^{\lambda}=\int\frac{d^{3}q}{(2\pi)^{3}}\epsilon^{\lambda}_{ij}(\hat{k})q^{i}q^{j}\bigg{(}2\Phi_{q}\Phi_{k-q}+(\mathcal{H}^{-1}\Phi^{\prime}_{q}+\Phi_{q})(\mathcal{H}^{-1}\Phi^{\prime}_{k-q}+\Phi_{k-q})\bigg{)}\,.$
(3.2)

The expression for the source term $S_{k}^{\lambda}(\tau)$ involves the
polarization tensor $\epsilon_{ij}^{\lambda}$. Furthermore, the Bardeen
potential $\Phi_{k}$ is related to the curvature perturbation through the
following relation

$\Phi_{k}=\frac{2}{3}\Gamma(k\tau)\mathcal{R}_{k},$ (3.3)

where $\Gamma(k\tau)$ is the transfer function. The power spectrum of
gravitational waves, corresponding to the solution of equation (3.1), is given
by

$\mathcal{P}_{h}(k,\tau)=2\int_{0}^{\infty}dt\int_{-1}^{1}dsF(s,t)\mathcal{P}_{\mathcal{R}}(kv)\mathcal{P}_{\mathcal{R}}(ku),$
(3.4)

where $u=\frac{s+t+1}{2}$ and $v=\frac{t-s+1}{2}$. For details of $F(s,t)$
see[74]. After calculating the power spectrum of gravitational waves, the
observed energy density of produced gravitational waves can be obtained as
follows

$\Omega_{GW}h_{0}^{2}=\Omega_{r}h_{0}^{2}\bigg{(}\frac{g_{*}}{g_{*,e}}\bigg{)}^{\frac{1}{3}}\Omega_{GW,e}(f),$
(3.5)

with $\Omega_{r}h_{0}^{2}\simeq 4.2\times 10^{-5}$ in which
$h_{0}=H_{0}/100\text{km}^{-1}$ and $\Omega_{r}$ is the current radiation
energy density. Also, $f=ck/(2\pi)$ and the factor
$\big{(}\frac{g_{*}}{g_{*,e}}\big{)}^{\frac{1}{3}}$ accounts for the change in
the number of relativistic degrees of freedom from the end of radiation
dominated era to the present time. In Figure 3 and 4 , the results of the
model are presented in comparison with the NANOGrav and future data. In these
figures, we have fixed the sharpness parameter to a specific value but
compared the results for different values of $c_{s}$. When comparing different
values of $c_{s}$, it is evident that for $c_{s}=1/4$ and $c_{s}=1/2$, the
model becomes qualitatively more consistent particularly with the first lines
of the NANOGrav data, as they cross the violin diagrams, which have a higher
probability according to the analysis in [1]. Conversely, $c_{s}=3/4$ lies
outside the violin diagrams for the first lines and does not explain the data
better than other values of $c_{s}$. Moreover, as seen in the left panel of
Fig. 3, various values of $c_{s}$ present a slightly different prediction for
future data, as $c_{s}=1/2$ and $c_{s}=3/4$ lie in a wider range of
predictions. However, in Fig. 4, $c_{s}=1$ lies in a narrower region of future
data.

Figure 3: Comparison of the model’s results with NANOGrav and future data for
$h=-6$. It can be seen that $c_{s}=1/4$ has a better fit, particularly with
the first violins, compared to other values of sound speed. Moreover,
$c_{s}=1/2$ and $c_{s}=3/4$ lie in a more wider range of future predictions.
Figure 4: Comparison of the model’s results with NANOGrav and future data for
$h=-12$. We observe that $c_{s}=1/4$ has a better fit, particularly with the
first violins, compared to other values of sound speed. Furthermore, other
values of $c_{s}$ can explain the future observations better than $c_{s}=1$.

##  4 Conclusion

The study of stochastic gravitational waves has become a valuable tool for
investigating the physics of both the early and late universe. In this
particular work, we have focused on comparing the predictions of a non-
attractor model for gravitational waves with the NANOGrav data from PTAs.

The non-attractor model considered in this study is characterized by three
distinct phases of inflation. The first and third phases correspond to slow-
roll inflation, which is a well-established paradigm in inflationary
cosmology. However, the second phase is unique as it exhibits non-attractor
behavior. This non-attractor phase introduces additional dynamics and features
into the inflationary process.

By comparing the predictions of the non-attractor model with the NANOGrav
data, we aim to assess the model’s viability and its ability to explain the
observed gravitational wave signals. This analysis provides valuable insights
into the physics of the early universe and enhances our understanding of the
inflationary dynamics and the generation of primordial gravitational waves. To
compare our results with the data from PTAs, we have also taken into account
observations of the PBH formation parameter, $f_{\text{PBH}}$. Notably, our
findings align well with the constraints imposed by the $f_{\text{PBH}}$ data.

In this study, we have explored a range of values for the sound speed
parameter, $c_{s}$, and the transition parameter, $h$. By considering these
various scenarios, we aim to comprehensively investigate the implications of
different parameter choices and their impact on the predicted outcomes. From
the plotted results, it is evident that models with $c_{s}=1/2$ and
$c_{s}=1/4$, exhibit better consistency with the observational data especially
in the first lines. However, it is important to note that these comparisons
are currently qualitative, and in order to provide a more precise assessment,
a quantitative approach will be pursued in future studies.

Moreover, we have considered the parameter space which have a good prediction
for the value of $f_{PBH}$ in the range of NANOGrav frequency. One may
consider the cases in which the power spectrum doesn’t grow enough. This may
make the model more consistent with PTAs data for other values of $c_{s}$. On
the hand one can study the effects of loop corrections in this study as well
and see how loop corrections can affect the predictions of model for the PTAs
data. The other possiblity that can be studied is to consider the case in
which the mode functions are not initially in a Bunch-Davis vacuum and see how
the predictions change with different choices of vacuum.

To summarize, this study emphasizes that models with lower sound speed
parameter values, $c_{s}$, exhibit more agreement with PTAs data. Nonetheless,
a meticulous and comprehensive data analysis is required to achieve a more
thorough assessment of the parameter space within the model and make a better
fit to the PTAs data. Another aspect worth investigating is the utilization of
a broader range of sound speed and sharp transition parameter values to derive
more dependable results for comparisons. However, this will require additional
computational power. We intend to conduct these analyses in future research
endeavors.

In addition to the aforementioned comparisons, we have extended the analytical
results presented in [56] to accommodate arbitrary values of the sound speed
parameter, $c_{s}$. Furthermore, we have computed the non-Gaussianity
parameter, $f_{\text{NL}}$ in this model as well. As expected, we find that
the non-Gaussianity parameter approaches zero as the transition becomes
milder. This aligns with the intuition that a smoother transition would lead
to a more Gaussian distribution of primordial fluctuations. By generalizing
the analytical results and investigating the non-Gaussianity parameter, we
gain a deeper understanding of the statistical properties of the primordial
perturbations in our model.

Acknowledgments: We would like to express our deep gratitude to Hassan
Firouzjahi for his valuable discussions and meticulous comments, which greatly
contributed to the development of this work. We are also thankful to Alireza
Talebian, Mohammad Hossein Namjoo and Sina Hooshangi for their insightful
discussions, which have enriched our understanding of the subject matter. K.
A. would like to extend special thanks to the University of SISSA and ICTP for
their generous hospitality during the conference ”Cosmology in Miramare,
2023,” where this paper was in progress.

## References

  * [1] G. Agazie et al. [NANOGrav], Astrophys. J. Lett. 951, no.1, L8 (2023)
  * [2] D. J. Reardon et al, Astrophys. J. Lett. 951, no.1, L6 (2023)
  * [3] J. Antoniadis et al, [arXiv:2306.16214 [astro-ph.HE]].
  * [4] H. Xu, S. Chen, Y. Guo, J. Jiang, B. Wang, J. Xu, Z. Xue, R. N. Caballero, J. Yuan and Y. Xu, et al. Res. Astron. Astrophys. 23, no.7, 075024 (2023)
  * [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).
  * [10] A. Vilenkin, Phys. Lett. B 107, 47-50 (1981).
  * [11] R. R. Caldwell and B. Allen, Phys. Rev. D 45, 3447-3468 (1992).
  * [12] A. Vilenkin, Phys. Rev. D 23, 852-857 (1981).
  * [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] K. N. Ananda, C. Clarkson and D. Wands, Phys. Rev. D 75, 123518 (2007).
  * [17] D. Baumann, P. J. Steinhardt, K. Takahashi and K. Ichiki, Phys. Rev. D 76, 084019 (2007).
  * [18] E. Bugaev and P. Klimai, Phys. Rev. D 81, 023517 (2010).
  * [19] H. Assadullahi and D. Wands, Phys. Rev. D 79, 083511 (2009).
  * [20] L. Alabidi, K. Kohri, M. Sasaki and Y. Sendouda, JCAP 09, 017 (2012).
  * [21] R. g. Cai, S. Pi and M. Sasaki, Phys. Rev. Lett. 122, no.20, 201101 (2019).
  * [22] S. Pi and M. Sasaki, JCAP 09, 037 (2020).
  * [23] S. Balaji, G. Domenech and J. Silk, JCAP 09, 016 (2022).
  * [24] A. Talebian, S. A. Hosseini Mansoori and H. Firouzjahi, Astrophys. J. 948, no.1, 48 (2023).
  * [25] G. Domènech, Universe 7, no.11, 398 (2021).
  * [26] B. Carr, F. Kuhnel and M. Sandstad, Phys. Rev. D 94, no.8, 083504 (2016).
  * [27] B. Carr and F. Kuhnel, Ann. Rev. Nucl. Part. Sci. 70, 355-394 (2020).
  * [28] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, Class. Quant. Grav. 35, no.6, 063001 (2018).
  * [29] O. Özsoy and G. Tasinato, Universe 9, no.5, 203 (2023).
  * [30] C. T. Byrnes and P. S. Cole, [arXiv:2112.05716 [astro-ph.CO]].
  * [31] P. Ivanov, P. Naselsky and I. Novikov, Phys. Rev. D 50, 7173-7178 (1994).
  * [32] J. Garcia-Bellido and E. Ruiz Morales, Phys. Dark Univ. 18, 47-54 (2017).
  * [33] M. Biagetti, G. Franciolini, A. Kehagias and A. Riotto, JCAP 07, 032 (2018).
  * [34] H. V. Ragavendra, P. Saha, L. Sriramkumar and J. Silk, Phys. Rev. D 103, no.8, 083510 (2021).
  * [35] H. Di and Y. Gong, JCAP 07, 007 (2018).
  * [36] J. Liu, Z. K. Guo and R. G. Cai, Phys. Rev. D 101, no.8, 083535 (2020)
  * [37] S. Hooshangi, A. Talebian, M. H. Namjoo and H. Firouzjahi, Phys. Rev. D 105, no.8, 083525 (2022).
  * [38] S. Hooshangi, M. H. Namjoo and M. Noorbala, JCAP 09, 023 (2023).
  * [39] A. Ghoshal, A. Moursy and Q. Shafi, [arXiv:2306.04002 [hep-ph]].
  * [40] W. H. Kinney, Phys. Rev. D 72, 023515 (2005).
  * [41] M. J. P. Morse and W. H. Kinney, Phys. Rev. D 97, no.12, 123519 (2018).
  * [42] W. C. Lin, M. J. P. Morse and W. H. Kinney, JCAP 09, 063 (2019).
  * [43] M. H. Namjoo, H. Firouzjahi and M. Sasaki, Europhys. Lett. 101, 39001 (2013).
  * [44] J. M. Maldacena, JHEP 0305, 013 (2003).
  * [45] P. Creminelli and M. Zaldarriaga, JCAP 10, 006 (2004).
  * [46] Y. F. Cai, X. Chen, M. H. Namjoo, M. Sasaki, D. G. Wang and Z. Wang, JCAP 05, 012 (2018).
  * [47] J. Martin, H. Motohashi and T. Suyama, Phys. Rev. D 87, no.2, 023514 (2013).
  * [48] X. Chen, H. Firouzjahi, M. H. Namjoo and M. Sasaki, EPL 102, no.5, 59001 (2013)
  * [49] X. Chen, H. Firouzjahi, E. Komatsu, M. H. Namjoo and M. Sasaki, JCAP 12, 039 (2013)
  * [50] X. Chen, M. x. Huang, S. Kachru and G. Shiu, JCAP 01, 002 (2007)
  * [51] M. Akhshik, H. Firouzjahi and S. Jazayeri, JCAP 12, 027 (2015).
  * [52] S. Mooij and G. A. Palma, JCAP 11, 025 (2015).
  * [53] R. Bravo, S. Mooij, G. A. Palma and B. Pradenas, JCAP 05, 024 (2018).
  * [54] B. Finelli, G. Goon, E. Pajer and L. Santoni, Phys. Rev. D 97, no.6, 063531 (2018).
  * [55] S. Pi and M. Sasaki, Phys. Rev. Lett. 131, no.1, 011002 (2023).
  * [56] H. Firouzjahi and A. Talebian, [arXiv:2307.03164 [gr-qc]].
  * [57] J. Kristiano and J. Yokoyama, [arXiv2211.03395 [hep-th]].
  * [58] J. Kristiano and J. Yokoyama, [arXiv:2303.00341 [hep-th]].
  * [59] A. Riotto, [arXiv:2301.00599 [astro-ph.CO]].
  * [60] A. Riotto, [arXiv:2303.01727 [astro-ph.CO]].
  * [61] S. Choudhury, S. Panda and M. Sami, Phys. Lett. B 845, 138123 (2023).
  * [62] S. Choudhury, S. Panda and M. Sami, [arXiv:2303.06066 [astro-ph.CO]].
  * [63] H. Firouzjahi, [arXiv:2303.12025 [astro-ph.CO]].
  * [64] H. Firouzjahi and A. Riotto, [arXiv2304.07801 [astro-ph.CO]].
  * [65] H. Motohashi and Y. Tada, JCAP 08, 069 (2023).
  * [66] G. Franciolini, A. Iovino, Junior., M. Taoso and A. Urbano, [arXiv2305.03491 [astro-ph.CO]].
  * [67] G. Tasinato, Phys. Rev. D 108, no.4, 043526 (2023).
  * [68] H. Firouzjahi, Phys. Rev. D 108, no.4, 043532 (2023).
  * [69] J. Fumagalli, [arXiv:2305.19263 [astro-ph.CO]].
  * [70] Y. Tada, T. Terada and J. Tokuda, [arXiv:2308.04732 [hep-th]].
  * [71] J. Garcia-Bellido, M. Peloso and C. Unal, JCAP 12, 031 (2016).
  * [72] S. Young and C. T. Byrnes, JCAP 08, 052 (2013).
  * [73] C. Cheung, A. L. Fitzpatrick, J. Kaplan and L. Senatore, JCAP 02, 021 (2008).
  * [74] K. Kohri and T. Terada, Phys. Rev. D 97, no.12, 123532 (2018).
  * [75] M. Akhshik, H. Firouzjahi and S. Jazayeri, JCAP 07, 048 (2015).
  * [76] M. H. Namjoo, H. Firouzjahi and M. Sasaki, JCAP 12, 018 (2012) doi:10.1088/1475-7516/2012/12/018 [arXiv:1207.3638 [hep-th]].