cubic interaction with two powers of $\sigma$’s, motivated by its potential of
detectability. Indeed, simpler single-exchange bispectra with a single power
of $\sigma$ in the cubic interaction generically lead to a small observable
signal. First, the “Lorentz-covariant” combination
$\propto\rho\times\partial_{\mu}\varphi\partial^{\mu}\varphi\sigma$ from the
unitary gauge operator $\rho\times\delta g^{00}\sigma$ has a strength fixed by
symmetries to be proportional to the quadratic mixing $\rho$; the latter being
treated in perturbation theory in most analytical calculations, this enforces
the corresponding bispectrum to be small. Even in the less explored regime of
a large quadratic mixing, see Refs. [28, 31, 170, 90], this channel can at
most lead to non-Gaussianities of order roughly unity [171]. Second, although
from the effective field theory point of view there may exist another cubic
vertex leading to a single-exchange bispectrum with a size independent from
the quadratic mixing (of the form
$\propto\tilde{\rho}\times(\varphi^{\prime})^{2}\sigma$ with $\tilde{\rho}$
potentially large), this interaction is not found in the large class of
general non-linear sigma models of inflation. On the contrary, as we have just
argued in the two previous paragraphs, cubic interactions leading to double-
exchange diagrams are both natural to consider from the effective field theory
point of view and found in practice in concrete multifield realizations of
inflation. The size of this interaction is not fixed by the small quadratic
mixing and may therefore be large while respecting the perturbative scheme
used for the derivation of our analytical results. And indeed, it was shown
that this channel could lead to potentially large non-Gaussianities both in
the small and large mixing regimes [171].

##### Perturbativity bound, naturalness and size of the bispectrum.

Although we have just argued that the signal from the double-exchange diagram
may be large, it cannot be arbitrary large. First, cubic interactions
naturally define strong coupling scales, ensuring that the perturbative
expansion in powers of fields’ fluctuations is a meaningful series. To be
precise, the cubic vertex responsible for the double-exchange diagram is
particular in the sense that $\lambda\varphi^{\prime}\sigma^{2}$ is a marginal
operator (of mass dimension 4 in 4 dimensions), so that the coupling $\lambda$
appearing in the time integrals is dimensionless and cannot be used to define
a strong coupling scale. However, for consistency of the perturbative
expansion in a similar fashion to the interpretation of strong coupling
scales, we should still ask that this cubic interaction remains smaller than
the kinetic (and mass) terms used to define the free theory [171], from which
one can derive the perturbativity bound

$\text{Perturbativity bound:}\quad|\lambda|\lesssim 1\,.$ (5.2)

Additionally, one may want to ask the size of the cubic interaction to be
natural in order to avoid requiring fine-tuning to explain the smallness of
loop corrections. In particular for the interaction under study, it was shown
in Ref. [171] that asking a mass of $\sigma$ naturally of order $H$ in the
effective field theory of multifield fluctuations framework, amounts to
setting 666 More precisely, this naturalness criterion corresponds to asking
that loop corrections to the mass of the heavy field, from such cubic
interactions, are not larger than the bare mass $m$ used in the calculation.
For $m\sim H$, this gives the upper bound (5.3), see Ref. [171]. Indeed,
naturalness is a statement about limiting the amount of fine-tuning required
to explain the value for a quantity given by the algebraic sum of several
terms. In the reasoning above, the fine-tuning for a mass of order Hubble with
a large $\lambda$ is clear: both the bare mass and the loop corrections should
be large and almost cancel to give a final value $\sim H$. Fine-tuning for a
large $\lambda$ can also be seen by specifying to concrete models of
inflation. For the large class of non-linear sigma models, the masses of
isocurvature fluctuations go as [174, 175] $m^{2}/H^{2}=-\lambda
f_{\pi}^{2}/H^{2}+V_{;ss}/H^{2}$, where $V_{;ss}$ represents consistent
projections of the Hessian matrix of the potential. Clearly, having a large
$\lambda$ while maintaining $m\sim H$ requires some degree of fine-tuning.
Because the corresponding naturalness criterion, $|\lambda|\lesssim
2\pi\Delta_{\zeta}$ for a weak quadratic mixing, is (slightly) more model-
dependent, we will not use it in the following.

$\text{Naturalness
criterion:}\quad|\lambda|\lesssim(2\pi\Delta_{\zeta})^{1/2}\,.$ (5.3)

Given the observed amplitude of curvature fluctuations at CMB scales, the
bound from naturalness is $\mathcal{O}(100)$ times more constraining than the
one from perturbativity. Note however that the status of the naturalness
criterion is less firm, as asking naturalness is not a prerequisite for a
theory to be well-behaved. Indeed, (approximate) cancellations may be
explained by (approximate) symmetries, accidental or not, or simply by an
appropriate amount of fine-tuning. On the contrary, perturbativity is not
negotiable in our framework, and we will always enforce it.

One can read the parametric dependence of the bispectrum for the inflaton’s
fluctuations $\varphi$, e.g. in Eq. (4.97), and translate it in terms of the
one for the curvature fluctuation. More precisely, we will be using the
dimensionless bispectrum shape function $S$, related to the bispectrum of
$\varphi$, in the weak mixing regime where $\Delta_{\varphi}=H/(2\pi)$, as
follows:

$\langle\varphi_{\bf{k}_{1}}\varphi_{\bf{k}_{2}}\varphi_{\bf{k}_{3}}\rangle^{\prime}=-\frac{H^{5}}{f_{\pi}^{2}(k_{1}k_{2}k_{3})^{2}}S(k_{1},k_{2},k_{3})\,.$
(5.4)

Without entering yet into the complicated mass and kinematical dependence of
the entire shape function, we read its parametric dependence related to the
size of the usual $f_{\mathrm{NL}}$’s parameters:

$S\sim
f_{\mathrm{NL}}\sim(\Delta_{\zeta})^{-1}\times(\rho/H)^{2}\times\lambda\,,$
(5.5)

where we used that $f_{\pi}^{2}$ is related to $\Delta_{\zeta}^{-1}$ in the
regime of validity of our analytical calculations. This confirms that, even in
the weak mixing regime with $\rho/H\ll 1$ and respecting the perturbativity
bound $|\lambda|\lesssim 1$, it is possible to have non-Gaussianities of order
one or more thanks to the large factor
$(\Delta_{\zeta})^{-1}\sim\mathcal{O}(10^{4})$. If one adds the requirement of
naturalness, it is still possible to get order one non-Gaussianities by
pushing towards intermediate values of the quadratic mixing,
$\rho/H\sim\mathcal{O}(10^{-1})$, or even beyond (a regime not encompassed by
our analytical approach from the bootstrap perspective but nevertheless
perfectly viable theoretically).

##### Single-field effective theory.

When the mass of the additional fluctuation $\sigma$ is large enough, it may
be integrated out of the theory [184, 185, 186], leading to an effectively
single-field theory for the curvature fluctuation only. This procedure
including all cubic interactions for two-fields non-linear sigma models was
performed in [174] and for any number of fields in [171]. Repeating the
procedure only for the cubic interaction considered here, we find the
following single-field effective theory in terms of the curvature fluctuation
directly:

$\displaystyle\mathcal{L}_{\mathrm{EFT}}$ $\displaystyle=\frac{a^{3}\epsilon
M_{\mathrm{Pl}}^{2}}{c_{s}^{2}}\left[\dot{\zeta}^{2}-c_{s}^{2}\frac{(\partial\zeta)^{2}}{a^{2}}+\left(\frac{1}{c_{s}^{2}}-1\right)A\frac{\dot{\zeta}^{3}}{H}\right]\,,\quad\text{with}$
(5.6) $\displaystyle\frac{1}{c_{s}^{2}}-1$
$\displaystyle=\left(\frac{\rho}{H}\right)^{2}\times\frac{H^{2}}{m^{2}}\quad\text{and}\quad
A=-2c_{s}^{2}\times\frac{\lambda
f_{\pi}^{2}}{H^{2}}\times\frac{H^{2}}{m^{2}}\,.$ (5.7)

Within this single-field theory, it is straightforward to compute the
bispectrum shape function in all kinematical regimes:

$S_{\mathrm{EFT}}(k_{1},k_{2},k_{3})=\frac{3A}{2}\left(\frac{1}{c_{s}^{2}}-1\right)\frac{k_{1}k_{2}k_{3}}{k_{123}^{3}}\,,$
(5.8)

where $k_{123}=k_{1}+k_{2}+k_{3}$. This shape not only peaks on equilateral
configurations but is also very correlated with the equilateral shape
template, consisting of:

$S_{\mathrm{eq}}(k_{1},k_{2},k_{3})=\left(\frac{k_{1}}{k_{2}}+\text{5
perms}\right)-\left(\frac{k_{3}^{2}}{k_{1}k_{2}}+\text{2 perms}\right)-2\,.$
(5.9)

Indeed, defining as usual the cosine between two shapes as their
correlation[187, 188]

$\displaystyle\mathrm{Cos}\left(S_{1},S_{2}\right)$
$\displaystyle=\frac{S_{1}\cdot S_{2}}{\sqrt{(S_{1}\cdot S_{1})(S_{2}\cdot
S_{2})}}\quad\text{with}$ (5.10) $\displaystyle S_{1}\cdot S_{2}$
$\displaystyle=\int_{1/2}^{1}\differential
x_{2}\int_{1-x_{2}}^{x_{2}}\differential x_{3}\,S_{1}(k,x_{2}\times
k,x_{3}\times k)S_{2}(k,x_{2}\times k,x_{3}\times k)\,,$ (5.11)

one finds $\mathrm{Cos}\left(S_{\mathrm{EFT}},S_{\mathrm{eq}}\right)\simeq
0.988$. We thus safely dub its amplitude at equilateral configurations the
parameter $f_{\mathrm{NL}}^{\mathrm{eq}}$ that one can look for in the data
with an equilateral shape template. To put it in a nutshell, the single-field
effective theory for the double-exchange diagram predicts an equilateral
bispectrum shape with amplitude

$f_{\mathrm{NL}}^{\mathrm{eq}}=-\frac{1}{9}\left(\frac{\rho}{H}\right)^{2}\frac{1}{1+\rho^{2}/m^{2}}\times\frac{\lambda
f_{\pi}^{2}}{H^{2}}\times\frac{H^{4}}{m^{4}}\,.$ (5.12)

Note that this single-field EFT does not rely on the assumption of a weak
quadratic mixing and is therefore valid for values of $\rho/H$ larger than
unity. In the regime of a weak quadratic mixing though (that we read actually
extends up to $\rho\lesssim m$), we recover the scaling with parameters
naively derived in Eq. (5.5), together with additional information in this
regime: the kinematical dependence of the shape is dominantly equilateral and
the exact prefactor including mass dependence is $-H^{4}/(18\pi m^{4})$. Note
however that the single-field effective theory relies on the hierarchy
$m^{2}c_{s}^{2}\gg H^{2}$ [174], and therefore implicitly also on the absence
of any other hierarchy in the theory. A subtlety comes about when looking at
soft limits of correlation functions, at which there is a kinematical
hierarchy. As far as the bispectrum is concerned, the unique soft limit is the
squeezed one. Let us fix $k_{3}$ the smallest momentum, in the squeezed limit
the ratio $k_{3}/k_{1}$ is very small and, therefore, one may well have
$(k_{3}/k_{1})\times(m^{2}c_{s}^{2}/H^{2})\lesssim 1$ although the field
$\sigma$ is very massive. Said otherwise, we expect the single-field
description to fail to capture not only the region of parameter space with
masses of order $H$ or smaller, but also the kinematical region corresponding
to soft limits. As well known, these soft limits indeed encode striking
signatures of particle production which cannot be described within the single-
field framework.

##### Exact numerical evolution.

The double-exchange bispectrum diagram may also be computed exactly from first
principles by computing its time evolution in the bulk of the inflationary
spacetime. It is precisely the aim of the so-called Cosmological Flow [170,
171] to propose a framework for automatically solving the differential
equations verified by the different two- and three-point correlators of any
inflationary theory defined at the level of fluctuations. The numerical
implementation of this approach is called CosmoFlow, it is open source and
available on GitHub and the paper describing it is Ref. [172]. It is
straightforward to implement the interactions leading to the double-exchange
bispectrum diagram within it, and in the coming sections we use it as a non-
trivial check of our analytical results. Indeed, importantly, the Cosmological
Flow approach does not rely on a perturbative scheme for the quadratic mixing.
Strictly speaking, it therefore does not solve for the same theory as the one
described in this work, however in the common regime of validity we should
expect both bootstrap analytical predictions and numerical calculations to
agree.

In the next sections, we explore the parametric and kinematical dependencies
of the double-exchange bispectrum. We then compare them to the ones of the
single-exchange channel.

###  5.1 Equilateral value and mass dependence

We start by exploring the size of the signal on the equilateral configuration
$k_{1}=k_{2}=k_{3}$ and its dependence on the mass of the double exchanged
field. At this particular point in kinematical space, we expect the single-
field effective description to be valid as soon as $m^{2}c_{s}^{2}/H^{2}\gg
1$. In practice we will focus on masses $m\gtrsim H$ and on a quadratic
coupling $\rho/H$ smaller than unity as required by consistency of the
bootstrap approach, a situation for which $c_{s}^{2}\simeq 1$.

In the left panel of Fig. 4, we compare the analytical prediction from the
multifield bootstrap approach with the single-field effective theory
description of Eq. (5.12) and the complete numerical resolution not relying on
a perturbative quadratic mixing with CosmoFlow, for $\rho/H=0.1$. The
agreement between the two exact methods in this regime is impressive. As an
example, both predict a change from positive to negative values with $\mu$
growing, precisely at the same threshold $\mu\simeq 1.5$. As expected, the
single-field effective theory prediction is only faithful at large values of
the mass, roughly $\mu\gtrsim 5$, and always predicts a negative bispectrum at
the equilateral configuration. In this regime of a small quadratic mixing and
a large mass, a simple fitting formula for our bootstrap result can be
derived,

$S_{\mathrm{bootstrap}}(k,k,k)\underset{\mu\gg
1}{\simeq}-0.11\times\left(\frac{\rho}{H}\right)^{2}\times\frac{\lambda
f_{\pi}^{2}}{H^{2}}\times\frac{1}{\mu^{4}}\,,$ (5.13)

which precisely matches the single-field effective theory result (5.12) with
the same numerical coefficient $-1/9\simeq-0.11$ (we remind that
$\mu=\sqrt{m^{2}/H^{2}-9/4}\simeq m/H$ in the regime of a large mass).

In the right panel of Fig. 4, we show the true size of the bispectrum at
equilateral configurations for $\rho/H=0.1$ as a function of the mass
parameter $\mu$ and the cubic coupling constant $\lambda$. We compare the
result to observational constraints from Planck,
$f_{\mathrm{NL}}^{\mathrm{eq}}=-26\pm 47$ at $1\sigma$ [189], by assuming that
the bispectrum can be qualitatively constrained with an equilateral shape, an
assumption to which we will come back in the next section. We find instructive
to see that the necessary theoretical constraint from perturbativity is close
to current observational constraints. We also show how, if one additionally
asks for naturalness, this reduces even more the possible size of the
bispectrum to less than unity.

(a)

(b)

Figure 4:  Left: Comparison of the values of the double-exchange bispectrum at
equilateral configurations from the exact numerical resolution using CosmoFlow
(blue line), the bootstrap analytical prediction of this work (orange dashed
line), as well for the corresponding single-field effective theory (green
line). The quadratic coupling is fixed as $\rho=0.1H$ and the mass of the
double exchanged field varies along with the horizontal axis. The time
evolution with CosmoFlow starts being numerically more demanding for larger
values of $\mu$ so it has been cut at $\mu\sim 10$ for purely computational
costs reasons. Note that the overall amplitude has been rescaled by the
constant factor $\lambda(2\pi\Delta_{\zeta})^{-1}$ that may be way larger than
unity, leading to potentially large primordial non-Gaussianities even in this
regime of a small quadratic mixing. Right: True size of the double-exchange
bispectrum at the equilateral configuration for $\rho/H=0.1$ as a function of
the mass parameter $\mu$ and the coupling constant $\lambda$. White regions
correspond to unobservably small values $<10^{-3}$ removed for clarity. The
contour line “observations” (solid line) stands for the $1\sigma$ constraint
on $f_{\mathrm{NL}}^{\mathrm{eq}}$ from Planck [189], while the other ones are
theoretical constraints. The perturbativity one (dashed line) is necessary to
enforce, while the naturalness one (dotted-dashed line) may or not be
required. The vertical valley visible by eye corresponds to the dip in the
left panel at which the bispectrum equilateral value changes from positive to
negative, from left to right in the figure, when $\lambda>0$.

###  5.2 Shape dependence

Our analytical expression for the result of the bootstrap equation
corresponding to the double-exchange bispectrum is valid for all kinematical
configurations of physical relevance. In particular, it is valid beyond the
equilateral configuration whose size dependence on the quadratic and cubic
couplings, as well as the mass of the exchange field, has been studied in the
previous section. Here we instead investigate the shape dependence on
$(x_{2},x_{3})$ of the bispectrum $S(k,x_{2}k,x_{3}k)$, normalized to unity on
equilateral configurations. We fix $\rho/H=0.1$ for definiteness, although in
the regime of weak mixing the shape is independent of it. In Fig. 5, we
compare the exact numerical prediction using CosmoFlow to the analytical
predictions derived from the bootstrap equation, for a few values of the mass
of the exchanged field. Unfortunately, the analytical formula converges slowly
in the folded limit region, corresponding to $x_{2}+x_{3}\simeq 1$ (the lower
ridge of the triangle). In practice we have therefore evaluated the shape on
the following kinematical region:
$x_{2}\in[0.58,1]\,,x_{3}\in[1.15-x_{2},x_{2}]$, which corresponds to roughly
$72\%$ of the total area. Note that the numerical resolution does not suffer
from this and can be evaluated quickly on any configuration but the very
squeezed ones corresponding to $x_{2}\sim 1$ and $x_{3}\ll 1$ (left corner of
the triangle), but for making the comparison explicit we evaluated it on the
same reduced triangle region as the analytical prediction. Agreement between
the two methods is impressive.

For completeness, we also calculate the correlations of the double-exchange
bispectrum shapes for a few representative values of the mass with the four
shape templates that are most used in data analysis: equilateral, flattened,
orthogonal and local bispectra. For this, we use the Cosine expression of Eq.
(5.10) and the shapes from CosmoFlow since we already showed agreement with
the bootstrap ones and that it allows to probe a larger kinematical region; we
simply cut at $x_{3}=0.1$ to avoid lengthy numerical evaluations in the
squeezed limit. The results are displayed in Fig. 6 below, where we also show
for reference the correlations between the shape templates for our kinematical
region. As expected, the shape resembles the local one for not-so-large masses
$m\lesssim H$ and the equilateral one for masses slightly larger than Hubble,
$m\gtrsim H$. The limiting case $m=3H/2$ seems to support almost the same
correlation with the local and the equilateral shape. However, we remind that
our correlations are computed in a restricted kinematical regime,
$x_{3}\geqslant 0.1$, and that having access to even more squeezed
configurations would reveal that all shapes are less correlated with the local
template. Even the case $m=H$ would be affected, as indeed the correct
template is the one of a quasi-local shape as first found in Ref. [5].

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: Comparison of the shape functions for three masses of order the
Hubble parameter obtained with CosmoFlow and the bootstrap method, for a
quadratic mixing $\rho/H=0.1$. The projected surface is shown for
representation purposes only and contains the same information as the
3D-shape.

Mass value | $m=H$ | $m=3H/2$ | $m=2H$ | $m=5H/2$ | $m=3H$ | $m=7H/2$
---|---|---|---|---|---|---
$\mathrm{Cos}(S,S_{\mathrm{eq}})$ | $0.71$ | $0.87$ | 0.87 | 0.92 | 0.93 | 0.93
$\mathrm{Cos}(S,S_{\mathrm{flat}})$ | $0.92$ | $0.84$ | $0.29$ | $0.72$ | $0.69$ | $0.72$
$\mathrm{Cos}(S,S_{\mathrm{orth}})$ | $-0.66$ | $-0.48$ | $0.14$ | $-0.32$ | $-0.28$ | $-0.31$
$\mathrm{Cos}(S,S_{\mathrm{loc}})$ | 0.99 | 0.90 | $0.36$ | $0.76$ | $0.75$ | $0.78$

$S_{i}$ | $S_{\mathrm{eq}}$ | $S_{\mathrm{flat}}$ | $S_{\mathrm{orth}}$ | $S_{\mathrm{loc}}$
---|---|---|---|---
$\mathrm{Cos}(S_{i},S_{\mathrm{eq}})$ | 1 | $0.46$ | $0.01$ | $0.66$
$\mathrm{Cos}(S_{i},S_{\mathrm{flat}})$ |  | 1 | $-0.88$ | $0.91$
$\mathrm{Cos}(S_{i},S_{\mathrm{orth}})$ |  |  | 1 | $-0.68$
$\mathrm{Cos}(S_{i},S_{\mathrm{loc}})$ |  |  |  | 1

Figure 6: Shape correlations between the double-exchange bispectrum at weak
mixing, $\rho/H=0.1$, with equilateral, flattened, orthogonal and local
bispectrum templates for a few representative values of the mass of the double
exchanged fluctuation. We also show for reference the correlations between
these shape templates for the kinematical region that we consider, with a cut
at $x_{3}=0.1$ in the squeezed limit. The dominant component is highlighted in
bold.

###  5.3 Cosmological collider signal

The CC signal lies in the soft limits of higher-order correlation functions,
chief amongst which is the squeezed limit in the primordial bispectrum. It
consists in an imprint left by the production of heavy particles during
inflation and is not encompassed by single-field effective theories. We have
already shown the explicit expression of the dominant squeezed bispectrum
corresponding to the exchange of two massive fluctuations of the same field
$\sigma$ in Eq. (4.97). Translated in terms of the primordial curvature
fluctuation bispectrum shape function, it reads:

$\displaystyle
S^{\mathrm{DE}}_{\mathrm{CC},\mathrm{LO}}=\left(\frac{\rho}{H}\right)^{2}\times\lambda(2\pi\Delta_{\zeta})^{-1}\times{\rm{Re}}\,\Bigg{\\{}\left[\frac{\pi^{1/2}}{2^{4+2\mathrm{i}\mu}}\frac{2\mathrm{i}\mu+5}{2\mu-3\mathrm{i}}\Gamma\left[\frac{1}{2}+\mathrm{i}\mu,-\mathrm{i}\mu\right]\big{(}1+{\rm{tanh}}(\pi\mu)\big{)}+\mathcal{O}\left(e^{-2\pi\mu}\right)\right]\kappa^{\frac{1}{2}+\mathrm{i}\mu}\Bigg{\\}}~{},$
(5.14)

where “DE” stands for Double-Exchange, “CC” for Cosmological Collider signal
as before, and “LO” for Leading-Order. We introduced the more usual squeezing
parameter $\kappa\equiv k_{1}/k_{3}$ and the squeezed limit corresponds to
$\kappa\ll 1$. Although this quantity (5.14) indeed represents the leading-
order squeezed behaviour, its mass suppression is exponential in the large
$\mu$ limit:
$S^{\mathrm{DE}}_{\mathrm{CC},\mathrm{LO}}\sim\mathcal{O}\left(\kappa^{1/2}\mu^{-1/2}e^{-\pi\mu}\right)$.
On the contrary the so-called background signal, denoted by “BG”, that is
decaying quicker in the squeezed limit, as $\kappa^{1}$, may be dominant in
not-so-squeezed configurations since its mass suppression is only polynomial:
$S^{\mathrm{DE}}_{\mathrm{BG}}\sim\mathcal{O}\left(\mu^{-4}\right)$ as can be
seen, for example, from its expression according to the single-field EFT and
that we derived in Eq. (5.12).

For practical predictions at intermediate squeezed configurations, we find it
useful to define and use next-to-leading-order (NLO) and next-to-next-to-
leading-order (NNLO) corrections to the LO cosmological collider signal shown
explicitly above. Let us clarify their meanings. In the final bootstrap
result, focusing on the squeezed limit, we are dealing with two different
series expansions. The first series is present in the analytical expressions
given by Eqs. (4.1)–(4.1)–(4.1), where we have a summation variable $m$. Each
term in the $m$ summation series is then expanded in terms of the squeezing
parameter $\kappa$. The terms “next-” and “next-to-next-” refer to the orders
in $\kappa$. For the CC signal part (containing non-analytic terms of the form
$\kappa^{\mathrm{i}\mu}$), the series expansion of the $m$-term in the sum
takes the following schematic form:
$S_{\mathrm{CC},m}^{\mathrm{DE}}\sim\kappa^{1/2+\mathrm{i}\mu+m}\left(1+\mathcal{O}(\kappa)+\mathcal{O}(\kappa^{2})+\ldots\right)$.
Moreover, explicit calculation shows that the $\mathcal{O}(\kappa)$ correction
to the LO $m=0$ term exactly cancels with the leading-order part of the $m=1$
term. Therefore, the total CC signal part scales as
$S_{\mathrm{CC}}^{\mathrm{DE}}\sim\kappa^{1/2+\mathrm{i}\mu}\left(1+\mathcal{O}(\kappa^{2})+\ldots\right)$.
Given the large suppression of the first non-zero correction to the CC signal
in the squeezed limit, it is clear that we should take into account the BG
part (containing only analytical terms, i.e. polynomial in $\kappa)$. However,
for the latter, each term in the $m$-summation begins at order $\kappa$ in the
squeezed limit:
$S_{\mathrm{BG},m}^{\mathrm{DE}}\sim\kappa^{1}\left(1+\mathcal{O}(\kappa)+\mathcal{O}(\kappa^{2})+\ldots\right)$.
In principle, to be consistent, one should therefore perform the sum over all
values of $m$. However, the higher the order in the $m$-summation, the more
they are suppressed by the mass parameter. Therefore, in the large mass limit,
$\mu\gtrsim 1$, it is consistent to keep only the first $m$-terms. In
practice, we therefore keep the first three terms corresponding to
$m=\\{0,1,2\\}$ to approximate the BG signal, and consistently keep the
dominant terms in the squeezed limit. The NLO correction corresponds to the
addition of the leading-order scaling $\mathcal{O}(\kappa)$ from this
approximated BG signal. The NNLO correction corresponds to the addition of the
$\mathcal{O}(\kappa^{2})$ terms from the BG signal. Note that the NNNLO one
would be the correction to the CC signal scaling as $\kappa^{5/2}$, but we do
not consider it explicitly here. Given the discussion above, it should be
clear that we expect the NLO and NNLO corrections to badly encapsulate the
cases where the mass is not large enough. Indeed, this can be seen for example
in the case $m=2H$ in Fig. 7, for which the LO term alone performs better. The
NLO and NNLO expressions are derived from our general bootstrap result, and
remain fully analytic. They are quite long and not enough illuminating to be
shown here, however they do allow for instantaneous evaluation for any value
of $\kappa$ and $\mu$, in opposition with the general result for which
squeezed configurations are computationally expensive to evaluate. As
mentioned previously, there exist potential methods to speed up the
convergence of the exact solution around folded/squeezed regions. For
instance, one approach involves employing the expansion series of the
hypergeometric function around the unit [177] to refine the results. However,
this transformation introduces considerable complexity, especially with
additional layers of series. Further transformation and simplification is
beyond the scope of this paper, and we leave this interesting discussion to
future work.

We now discuss the relative size of the cosmological collider signal from the
double-exchange channel that we have been focusing on in this work, compared
to the one from the single-exchange channel. For concrete comparison, we
consider the only cubic interaction leading to a single-exchange diagram that
is encountered in non-linear sigma models of inflation. We have already
mentioned that its strength is fixed by a non-linearly realized symmetry to be
proportional to the quadratic mixing,
$\rho/(2f_{\pi}^{2})\times\partial_{\mu}\varphi\partial^{\mu}\varphi\sigma$
(see Eq. (2.45)). The explicit expression of the cosmological collider signal
from this single-exchange interaction, first found in Ref. [4] and then
exactly confirmed with an independent method in Ref. [77], is shown in Eq.
(5.18) below. We find the ratio of the two squeezed limits, at leading-order
in $\kappa$ and $\mu$, to be remarkably simple:

$\frac{S^{\mathrm{DE}}_{\mathrm{CC},\mathrm{LO}}}{S^{\mathrm{SE}}_{\mathrm{CC},\mathrm{LO}}}=\lambda(2\pi\Delta_{\zeta})^{-1}\times\frac{2}{\mu^{2}}+\mathcal{O}\left(\mu^{-3}\right)\,,$
(5.15)

where “SE” denotes this Single-Exchange channel. The parametric suppression in
terms of the mass is exactly the same as the relative size of the equilateral
background signals from both channels (indeed
$S^{\mathrm{SE}}_{\mathrm{CC},\mathrm{LO}}\sim\mathcal{O}\left(\kappa^{1/2}\mu^{3/2}e^{-\pi\mu}\right)$
and $S^{\mathrm{SE}}_{\mathrm{BG}}\sim\mathcal{O}\left(\mu^{-2}\right)$). In
particular, the absence of further “Boltzmann” exponential suppression was
guessed in [125], and here we have rigorously proved it. Importantly, this
implies that for values of the mass leading to a signal potentially
observable, corresponding to $\mu$ of at most a few, the cosmological collider
signal from the double exchange bispectrum is dominant as long as
$\lambda(2\pi\Delta_{\zeta})^{-1}\gtrsim 1$. Such values are perfectly allowed
by perturbativity bounds, and can even be natural, as discussed at the
beginning of Sec. 5, therefore making the double-exchange channel an
interesting prospect for detecting the cosmological collider signal.

We now compare our exact bootstrap results (limited to close-to-equilateral
configurations $\kappa\gtrsim 0.1$) with the LO, NLO and NNLO expansions in
the squeezed limit, as well as the exact numerical results with CosmoFlow
(whose only limitation is the computational cost of very squeezed
configurations, in practice we restrict to $\kappa\gtrsim 10^{-3}$). The
results are shown in Fig. 7, in the weak mixing regime for $\rho/H=0.1$ where
our analytical predictions hold, and for a few representative values of the
mass of the massive field. It is interesting to notice the excellent agreement
between CosmoFlow results and i) the exact bootstrap result where it can be
evaluated efficiently, close to equilateral configurations; ii) the LO
expansion in the squeezed limit; iii) the NLO and NNLO expansions for not-so-
squeezed configurations when the mass is sizeable ($\mu\gtrsim 1)$.
Altogether, our example showcases the utility and synergy of analytical and
numerical methods, with precise and trustworthy predictions for the double-
exchange bispectrum for all configuration values, ranging from equilateral to
extremely squeezed ones.

(a)

(b)

(c)

(d)

Figure 7:  Cosmological collider signals in isosceles triangle configurations
of the primordial bispectrum, for the double-exchange channel in the weak
mixing regime with $\rho/H=0.1$, and normalized to unity in equilateral
configurations, for four representative values of the mass $m$.

For completeness, we will also quickly explore the strong quadratic mixing
regime. The phenomenology of the cosmological collider signal in this regime
has been thoroughly explored numerically in [171]. One of the main lessons is
the rescaling of the mass to an effective mass by the large quadratic mixing:
$m^{2}\rightarrow m_{\mathrm{eff}}^{2}=m^{2}+\rho^{2}$, visible in the
frequency of the oscillations in the cosmological collider signal, now
dictated by $\mu_{\mathrm{eff}}=\sqrt{\mu^{2}+\rho^{2}/H^{2}}$. In Fig. 8, we
explore how the analytical formula performs under such substitution compared
to the exact numerical result obtained in the strong mixing regime for
$\rho/H=1$ with CosmoFlow, though the analytical formulas are valid only for
the weak quadratic mixing regime by construction. Note that for this value of
the quadratic mixing, the true size of the bispectrum is exactly 100 times
bigger than the one shown in the right panel of Fig. 4, pushing it to a level
comparable to current and upcoming observational constraints, even for natural
values of the cubic coupling constant $\lambda$ that are well below the
perturbativity bound, thus making an impressive case for the double-exchange
channel. Interestingly, the modified analytical formulas correctly capture not
only the frequency of the oscillations, but also the overall shape of the
signal. They seem to miss however a slight damping of the CC signal, as well
as a small phase shift towards larger values. It would be interesting to
pursue analytical calculations tailored to the strong quadratic mixing regime
and precisely predict these interesting features, but this goes beyond the
scope of this article.

(a)

(b)

Figure 8:  Cosmological collider signals in isosceles triangle configurations
of the primordial bispectrum, for the double-exchange channel in the strong
mixing regime with $\rho/H=1$, and normalized to unity in equilateral
configurations, for two representative values of the effective mass
$m_{\mathrm{eff}}$. The bootstrap predictions are obtained by the naive
replacement $\mu\rightarrow\mu_{\mathrm{eff}}$ in our analytical formulas.
This replacement works surprisingly well, even though it misses a slight
damping of the CC signal as well as a small phase shift towards larger values.

###  5.4 Disentangling double-exchange from single-exchange channels

We have shown that the bispectrum signal from the double-exchange channel can
easily be larger than the one from the single-exchange one. We proved this to
be correct both for the equilateral configuration value and for the
cosmological collider signal lying in the squeezed limit. However, we still
lack a striking feature of the double-exchange channel that cannot be mimicked
by the single-exchange one. In this section, we explore three directions in
order to address this question: i) information in the phase of the
cosmological collider signal; ii) the exchange of two massive fluctuations
pertaining to two different isocurvature species with unequal masses; iii) the
information in the primordial trispectrum.

##### Phase information in the CC signal.

The shape functions of single-exchange and double-exchange bispectra, on
isosceles squeezed configurations, can be written as:

$\displaystyle S_{\mathrm{CC},\mathrm{LO}}^{\mathrm{SE}}$
$\displaystyle=\left(\frac{\rho}{H}\right)^{2}\times\mathrm{Re}\left[\kappa^{1/2+\mathrm{i}\mu}\mathcal{A}_{\mathrm{SE}}(\mu)e^{\mathrm{i}\,\delta_{\mathrm{SE}}(\mu)}\right]\,,$
(5.16) $\displaystyle S_{\mathrm{CC},\mathrm{LO}}^{\mathrm{DE}}$
$\displaystyle=\left(\frac{\rho}{H}\right)^{2}\times\lambda(2\pi\Delta_{\zeta})^{-1}\times\mathrm{Re}\left[\kappa^{1/2+\mathrm{i}\mu}\mathcal{A}_{\mathrm{DE}}(\mu)e^{\mathrm{i}\,\delta_{\mathrm{DE}}(\mu)}\right]\,.$
(5.17)

Here $\mathcal{A}_{\mathrm{SE},\mathrm{DE}}$ denote the amplitudes of the CC
signals and $\delta_{\mathrm{SE},\mathrm{DE}}$ their phases. Importantly, in
this squeezed limit and in the regime of validity (weak quadratic mixing) of
the analytical calculations, they only depend on $\mu$, the mass parameter of
the unique isocurvature species. The expression for
$\mathcal{A}_{\mathrm{DE}}(\mu)e^{\mathrm{i}\,\delta_{\mathrm{DE}}(\mu)}$ can
be directly read off Eq. (5.14), while the one for the single-exchange channel
can be found explicitly in Refs. [4, 77]:

$\mathcal{A}_{\mathrm{SE}}(\mu)e^{\mathrm{i}\,\delta_{\mathrm{SE}}(\mu)}=-\mathrm{i}\frac{\pi^{3/2}}{2^{4+2i\mu}}\frac{\Gamma(7/2+\mathrm{i}\mu)}{(2\mu-\mathrm{i})\Gamma(1+\mathrm{i}\mu)}\tanh(\pi\mu)e^{-\pi\mu}\left[\coth(\pi\mu)-\mathrm{i}\csch(\pi\mu)+1\right]^{2}\,.$
(5.18)

In the regime of validity of this analysis, the frequency of the CC signal is
directly given by $\mu$. It was argued in Ref. [82] that the phase of the CC
signal may bring additional information to help distinguishing the channel at
the origin of this signal. In the latter reference, differences in the phases
were investigated between two cubic interactions leading to single-exchange
diagrams: the standard Lorentz-covariant interaction denoted as “SE” in our
work and a more exotic Lorentz-violating one with schematic form
$\tilde{\rho}\varphi^{\prime 2}\sigma$, but nevertheless allowed from an
effective field theory point of view. Here, we compare the phases between the
ubiquitous single-exchange “SE” channel with Lorentz-covariant interaction and
the double-exchange one for the first time. As can be seen in the left panel
of Fig. 9, information about the phase, once the frequency $\mu$ is
determined, indeed places interesting constraints on the possible channel to
generate the CC signal. At large values of $\mu$ though, the two phases become
indistinguishable, but fortunately this is any way the regime for which the CC
signal is unobservably small. In the right panel of Fig. 9, we instead
investigate the dependence of an effective phase
$\delta_{\mathrm{DE}}(\mu,\kappa)$ in mildly squeezed configurations, as a
function of $\kappa$ and for a few values of the mass, when adding the first
non-vanishing correction to the CC signal in the squeezed limit expansion,
scaling as $\kappa^{5/2}$ as already mentioned. These not-so-squeezed
configurations are of particular relevance to observations, and it is
therefore important to dispose of full predictions for all kinematic regimes
as we do.

Figure 9:  Left: Comparison of the phases of single- and double-exchange
diagrams, as defined in Eqs. (5.16) and (5.17). The figure shows the one-to-
one correspondences between frequency and phase of the oscillations, unveiling
the potentiality to disentangle different CC channels from the phase
information. This is particularly true whenever $\mu$ is not too large, which
is the most interesting case. Both dashed lines represent the leading
contribution in the large $\mu$ limit, with a single power of the Boltzmann
suppression, $\mathcal{O}(e^{-\pi\mu})$. Solid lines represent the total
result including all contributions, making for a non-negligible correction in
the small $\mu$ region. Right: Evolution of the phase, for three
representative masses, as a function of the squeezing parameter when adding
the first non-vanishing correction scaling as $\kappa^{5/2}$ in the CC signal.
Only the region $\kappa\in[0.1,1]$ is affected by the deviation from a
constant, and we expect higher-order corrections in $\kappa$ in the CC signal
to also play a role when $\kappa\rightarrow 1$.

##### Multiple isocurvature species: inflationary flavor oscillations.

Before the current paragraph, all phenomenological consequences that we have
shown were derived under the assumption that there exists an unique massive
species $\sigma$, with mass $m$ and mass parameter $\mu$. This choice was
motivated by the long-standing tradition to explore CC signals in this
restricted class of multiple fields content, and under the implicit assumption
that only the lightest of massive degrees of freedom may leave an imprint on
cosmological observables. However, it was discovered in Ref. [77] that even
heavier species may strongly impact the primordial bispectrum in its squeezed
limit, due to an interesting interplay between the values of the masses
($m_{\alpha})$ and of the mixing angles of the theory (in $\rho_{\alpha}$)
that define mass eigenstates with respect to the flavor ones. The most
striking manifestation of multiple massive species is the appearance of
modulated oscillations in place of the usual single-harmonics CC signal,
signaling the presence of inflationary flavor oscillations analogous to those
appearing in the Standard Model of particle physics. By determining the
frequencies and relative sizes of the different modes, it was proposed that
observations could help measuring not only the masses but also the mixing
angles of the inflationary theory. Explicitly, the sum of all single-exchange
diagrams with multiple isocurvature species gives the following leading-order
CC signal in the squeezed limit [77]:

$\displaystyle
S_{\text{CC,LO}}^{\text{SE,multi}}=\sum_{\alpha}\left(\frac{\rho_{\alpha}}{H}\right)^{2}\real\left[\kappa^{1/2+\mathrm{i}\mu_{\alpha}}\mathcal{A}_{\text{SE}}(\mu_{\alpha})e^{\mathrm{i}\delta_{\text{SE}}({\mu_{\alpha}})}\right]\,,$
(5.19)

where $\mathcal{A}_{\text{SE}}$ and $\delta_{\text{SE}}$ were already defined
in Eqs. (5.16) and (5.18). Here, we extend the discussion of the latter
reference to the double-exchange channel. Where previously we had the exchange
of two fluctuations of the same species $\sigma$, we now allow for multiple
species $\sigma_{\alpha}$ with a priori unequal mass parameters
$\mu_{\alpha}$.

In the case of multiple isocurvature species, the double-exchange diagram
results in a superposition of oscillatory signals in the squeezed limit, with
the dominant contribution being given by

$\displaystyle
S_{\text{CC,LO}}^{\mathrm{DE},\mathrm{multi}}=\sum_{\alpha,\beta}\frac{\rho_{\alpha}\rho_{\beta}}{H^{2}}\times\lambda_{\alpha\beta}(2\pi\Delta_{\zeta})^{-1}\times\real\quantity[\kappa^{1/2+\mathrm{i}\mu_{\alpha}}{\mathcal{A}_{\mathrm{DE}}^{\mathrm{multi}}(\mu_{\alpha},\mu_{\beta})}e^{\mathrm{i}\delta_{\mathrm{DE}}^{\mathrm{multi}}({\mu_{\alpha},\mu_{\beta}})}],$
(5.20)

where $\mathcal{A}_{\mathrm{DE}}^{\mathrm{multi}}$ and
$\delta_{\mathrm{DE}}^{\mathrm{multi}}$ are uniquely determined from our
analytical formulae (4.104). Disentangling individual contributions $\beta$ to
the same frequency signal set by $\mu_{\alpha}$ is impossible, so we also
define the effective amplitude and phase of individual modes in the CC signal
as:

$\displaystyle S_{\text{CC,LO}}^{\text{DE,multi}}$
$\displaystyle=\sum_{\alpha}\left(\frac{\rho_{\alpha}}{H}\right)^{2}\real\left[\kappa^{1/2+\mathrm{i}\mu_{\alpha}}\tilde{\mathcal{A}}_{\text{DE}}^{\alpha}e^{\mathrm{i}\tilde{\delta}_{\text{DE}}^{\alpha}}\right]\,,$
(5.21)
$\displaystyle\text{with}\quad\tilde{\mathcal{A}}_{\text{DE}}^{\alpha}e^{\mathrm{i}\tilde{\delta}_{\text{DE}}^{\alpha}}$
$\displaystyle=\sum_{\beta}\frac{\rho_{\beta}}{\rho_{\alpha}}\times\lambda_{\alpha\beta}(2\pi\Delta_{\zeta})^{-1}\times\mathcal{A}_{\mathrm{DE}}^{\mathrm{multi}}(\mu_{\alpha},\mu_{\beta})e^{\mathrm{i}\delta_{\mathrm{DE}}^{\mathrm{multi}}({\mu_{\alpha},\mu_{\beta}})}\,,$

where we have factorised out two factors of $\rho_{\alpha}$ for simpler
comparison to the single-exchange diagram. Although a detailed study of the
phenomenology of the multiple isocurvature species double-exchange
diagrams—and its comparison to the single-exchange one—is beyond the scope of
this article, we already note an interesting feature. When there is a single
isocurvature species, we have seen in Fig. 9 that the phase is uniquely fixed
by the frequency of the CC signal, for both kinds of diagrams. With several
isocurvature species, a similar statement still holds for the single-exchange
diagram: each individual mode $\alpha$ in the CC signal has a phase uniquely
fixed by its frequency. However, this one-to-one correspondance is broken in
the double-exchange case, as both the amplitude and the phase of the mode
$\alpha$ depend also on all other isocurvature species’s masses $\mu_{\beta}$
and mixing angles in $\rho_{\beta}$, as long as their cubic interactions as
encoded in $\lambda_{\alpha\beta}$ have non-trivial non-diagonal elements.
Therefore, in the occasion that we can measure both the effective amplitude
and phase of different modes in the CC signal, we could easily rule out
single-exchange diagrams if the consistency condition between phase and
frequency is not verified. The opposite is harder, as a single mode in a
double-exchange diagram signal can mimic the phase of a single-exchange one.
In Fig. 10, we illustrate these new effects from the non-diagonal elements of
the cubic interactions in $\lambda_{\alpha\beta}$. We plot the signal in the
squeezed limit corresponding to inflationary flavor oscillations for a model
with $\varphi$ coupled to two massive scalar fluctuations
$(\sigma_{1},\sigma_{2})$, with mass parameters $(\mu_{1},\mu_{2})=(1.5,0.5)$.
The mixing angle $\theta_{12}$ between the inflationary flavor and mass
eigenstates is defined by
$(\rho_{1},\rho_{2})=\left(\rho\cos(\theta_{12})\,,\rho\sin(\theta_{12})\right)$.

Figure 10: Inflationary flavor oscillations in the squeezed limit of the
bispectrum from the double-exchange diagram. In both the left and right
panels, the masses of different isocurvature species are set to $\mu_{1}=1.5$
and $\mu_{2}=0.5$, while the mixing angles $\theta_{12}$ are chosen as
$\pi/10$ (left panel) and $\pi/4$ (right panel). The four distinct cases,
corresponding to different forms of the cubic interaction matrix
$\lambda_{\alpha\beta}$, are described in the main text below this figure. For
better visualization, we have introduced a $\kappa^{-1/2}$ rescaling factor,
and the shape function is arbitrarily normalized to its value at
$\kappa=10^{-1}$ .

It was shown in Ref. [77] that this angle is found from the rotation matrix
used to diagonalize the non-trivial mass matrix present in the flavor basis.
The left panel corresponds to $\theta_{12}=\pi/10$, which was identified in
the previous work as leading to striking signatures even from the single-
exchange channel. The right panel corresponds to $\theta_{12}=\pi/4$, for
which the mode mixing in the single-exchange channel can hardly be seen. In
each of the panels, we vary the form of the cubic interaction matrix
$\lambda_{\alpha\beta}\equiv\lambda\times e_{\alpha\beta}$ as follows, where
we only mention the non-zero components: Case 1 has $e_{11}=1$ only and
forbids the presence of mode mixing; Case 2 has $e_{11}=e_{22}=1$ and
corresponds to a situation analogous to the single-exchange diagram (mode
mixing can only be seen for small mixing angles); Case 3 has
$e_{11}=e_{12}=e_{21}=1$ and Case 4 has $e_{11}=e_{22}=e_{12}=e_{21}=1$, they
both correspond to new features of the double-exchange channel not
reproducible by the single-exchange one. Although definitely leading to
signals of different natures in the squeezed limit, which can help to
disentangle single-exchange diagrams from double-exchange ones from
observations in principle, it would be desirable to identify another striking
signature of the exchange of two different isocurvature species. In this
perspective, we now turn to a first study of the CC signal in the trispectrum
from a double-exchange diagram.

##### Towards the primordial trispectrum.

In the trispectrum, as evidenced by the patterns observed in Figure 3, the
oscillation features are totally different from those of the single-exchange,
appearing as if two signals with different frequencies were superimposed.
Therefore, at the trispectrum level, in principle, we can clearly distinguish
between single-exchange and double-exchange through the oscillation patterns.
In this subsection, we will delve into the details of these distinct features.
Firstly, let us define the dimensionless trispectrum $\mathcal{T}$ as follows

$\displaystyle\langle\varphi_{{\bf{k}}_{1}}\varphi_{{\bf{k}}_{2}}\varphi_{{\bf{k}}_{3}}\varphi_{{\bf{k}}_{4}}\rangle^{\prime}=\frac{H^{8}}{f_{\pi}^{4}}\frac{(k_{1234}/4)^{3}}{(k_{1}k_{2}k_{3}k_{4})^{3}}\mathcal{T}({\bf{k}}_{1},{\bf{k}}_{2},{\bf{k}}_{3},{\bf{k}}_{4})~{}.$
(5.22)

The interactions chosen here are specified in (2.11) and (2.16), then the
four-point functions can be related to the seed integral through the simple
relation (2.18). With these coupling choices, the $\tau_{\text{NL}}$ typically
denoting the size of the trispectrum, has the following parametric dependence:

$\displaystyle\mathcal{T}\sim\tau_{\text{NL}}\sim\left({\rho}/{H}\right)^{2}\times\left(H^{2}\tilde{\lambda}\right)\times(2\pi\Delta_{\zeta})^{-2}~{}.$
(5.23)

The kinematic dependence of the trispectrum from this diagram is relatively
simple, with no dependence on the diagonals of the quadrangle formed in
Fourier space by the four wavevectors $\bf{k}_{i}$, so we will simply write
$\mathcal{T}(k_{1},k_{2},k_{3},k_{4})$. Similar to what we did in the
bispectrum case, in the regime of perturbative quadratic mixing and the large
mass, we can easily find a simple fitting formula from the analytical result
in regular configurations:

$\mathcal{T}(k,k,k,k)\underset{\mu\gg
1}{\simeq}0.07\times\rho^{2}\tilde{\lambda}\times\frac{f_{\pi}^{4}}{H^{4}}\times\frac{1}{\mu^{4}}~{}.$
(5.24)

To investigate the phenomenon indicated in Figure 3, we will focus on one
particular channel of the trispectrum, setting $k_{2}=k_{4}=k$ and
$k_{1}=k_{3}=rk$, with $r\rightarrow 0$ representing the double soft limit. To
gain a clear understanding of the origin of this phenomenon, instead of
dealing with the complicated exact solution, we will rely on some approximate
formulae discussed in Section 4.2. Specifically, the approximation here arises
from a series expansion around $r\sim 0$ of the exact solutions. Under this
limit, we divide the approximate expression into three distinct parts, each
with their own interpretations,

$\displaystyle\mathcal{T}_{\text{approx}}=\mathcal{T}_{\text{ss}}+\mathcal{T}_{\text{sb}}+\mathcal{T}_{\text{bb}}~{}.$
(5.25)

$\mathcal{T}_{\text{ss}}$ arises when both momenta associated with massive
propagators are soft, leading to a CC signal from each. Given that the time-
ordered part $\mathcal{I}_{+++}$ dominates due to lighter Boltzmann
suppression, let us simplify and focus solely on this term for clarity. In
this context, $\mathcal{T}_{\text{ss}}$ originates from the homogeneous
solution $\mathcal{Y}^{\text{-2}~{}0~{}\text{-2}}_{+++}$, under the double
soft limit, analogous to Eq. (4.76) discussed earlier. For
$\mathcal{T}_{\text{sb}}$, the CC signal is generated only from one
propagator, while the other contribution arises from the particular solution,
similar to Eq. (4.82) under the double-soft limit. Additionally, due to
permutations, there is another contribution where only one massive propagator
becomes soft, approximated by Eq. (4.2). The last term
$\mathcal{T}_{\text{bb}}$, represents the leading-order analytical background.

Figure 11:  The comparison between the numerical calculation and analytical result. The gray dots, blue solid line and red line represents numerical results, approximation expression and exact result from (4.1) respectively. In the figure we set $k_{2}=k_{4}=k$ amd $k_{1}=k_{3}=rk$, the masses are chosen as $\mu_{\alpha}=\mu_{\beta}=2$ and the coupling $\tilde{\lambda}\cdot\rho^{2}$ is set to $10^{-2}$. The left region (light red) is dominated by the term $r^{2\mathrm{i}\mu}$ while the right region (light blue) is dominated by the oscillation $r^{\mathrm{i}\mu}$. For better visualization, we have multiplied the dimentionless trispectrum $\mathcal{T}(rk,k,rk,k)$ by a factor $r^{-7/2}$.  different terms | mass dependence |  $r$ dependence
---|---|---
$\mathcal{T}_{\text{ss}}$ | $\mu^{3/2}e^{-2\pi\mu}$ | $r^{3+2\mathrm{i}\mu}$
$\mathcal{T}_{\text{sb}}$ | $\mu^{3/2}e^{-\pi\mu}$ | $r^{7/2+\mathrm{i}\mu}$
$\mathcal{T}_{\text{bb}}$ | $\mu^{-4}$ | $r^{4}$
Table 1: The leading order mass and $r$ dependence of each term in Eq. (5.25).

Because the formula is lengthy even for the leading-order series expansion, we
do not intend to present details of in this section, but it can be directly
obtained using the methods introduced in Section 4.2. Here, we summarize the
key information about the magnitude and dependence on $r$ of each term in
Table 1, where we have only kept the leading order in the series expansion
around $r=0$. Now, it is clear to see the origins of the unique oscillation
pattern. $\mathcal{T}_{\text{ss}}$ excites CC signals in both massive
propagators, thus the signal is suppressed by $\mathcal{O}(e^{-2\pi\mu})$, and
the frequency is doubled to $2\mu$. As a compensation, it features a dominant
scaling with respect to the soft variable $r$ when it approaches zero, and the
overall $r$ dependence is $r^{3+2\mathrm{i}\mu}$. In contrast,
$\mathcal{T}_{\text{sb}}$ is suppressed by $\mathcal{O}(e^{-\pi\mu})$ and
oscillates as $r^{7/2+\mathrm{i}\mu}$. Consequently, there is a competition
between mass and $r$ dependence. In very soft regions, the term with the lower
power in $r$ will dominate, whereas in intermediate regions, the term with
less Boltzmann suppression will prevail. To illustrate this idea, in Figure
11, we compare the analytical expression with the numerical results. The
approximation depicted by the blue line closely matches the numerical
calculations represented by the gray dots, except in the region close to the
equilateral configuration. In this region, we employed the exact solution
shown as the red line, which actually matches perfectly the numerical results.
In this figure, we have assumed equal masses with $\mu=2$ and we have
multiplied the dimensionless trispectrum $\mathcal{T}(rk,k,rk,k)$ by a factor
of $r^{-7/2}$ for better visualization. Now, based on our discussion above, we
can understand the superimposed oscillation shape in Figure 11. This shape
arises from the competition between the oscillation signals
$\mathcal{O}(e^{-2\pi\mu})\times r^{2\mathrm{i}\mu+3}$ and
$\mathcal{O}(e^{-\pi\mu})\times r^{\mathrm{i}\mu+7/2}$. In the very soft
region (with a light red background in Figure 11), the term from the
homogeneous solution dominates due to the additional $r^{-1/2}$ factor, and
the frequency is $2\mu$. In the blue regions, due to lighter Boltzmann
suppression, $r^{\mathrm{i}\mu+7/2}$ dominates with a frequency of $\mu$. This
explains why the frequency on the left side of Figure 11 is double that of the
right region, making for a definitely striking feature of the double-exchange
channel.

##  6 Conclusions and Outlooks

Cosmological correlators inherit information about all primordial fields, with
any masses, spins, mixing angles, and including all kinds of interactions with
the curvature perturbation and the two polarisations of the gravitational
waves. In order to have at hand full predictions in any given inflationary
scenario, as well as for not biasing the interpretation of (upcoming)
cosmological data, channels involving multiple massive exchanges must be taken
into account. However, obtaining analytical solutions for these channels is
highly challenging. In this work, we present the first example wherein we
apply recently developed bootstrap equations to find the exact and explicit
solutions of four-point and three-point correlation functions with the
exchange of two massive fields. The approach employed in this paper involves
applying differential operators associated with each massive propagator
individually, thereby expressing derivatives of the double-exchange seed
integral as the corresponding single-exchange one without this propagator. By
utilizing known results from prior research as input, we ultimately derive
exact analytical solutions for these bootstrap equations. Through proper
selection of variables, our final results exhibit a simplified form,
comprising only one layer of series summation. This simplification facilitates
faster convergence and makes it easier to investigate the analytic structure.
The primary results of the four-point function exhibit divergences in certain
terms for specific kinematic regions, making the task of taking the three-
point function limit challenging. However, we observed that by employing the
continuation of certain special functions and regrouping the divergences in
pairs, treating all contributions as a whole, the final expression is
explicitly convergent. We also conducted a detailed analysis of the
phenomenology related to double-exchange correlation functions, together with
a careful comparison with an independent numerical method using CosmoFlow. We
concluded that the primordial bispectrum generated from this channel may be
large, both at equilateral configurations and in the squeezed limit where the
cosmological collider dominates, thereby opening the thrilling possibility to
probe the existence of massive primordial fields via cosmological
observations. We explored diverse phenomenological aspects that make the
double-exchange channels particularly intriguing, including specific phase
information in the CC signal, new features in the inflationary flavor
oscillations pattern and a unique transition to double frequency in the double
soft limit of the primordial trispectrum. More generally, our work showcases
the utility of using diverse methods for the calculations of cosmological
correlators, and taking full advantage of their synergies.

There are many aspects that deserve further investigation in the future. Here
we outline a few of them. First, as previously mentioned in the main text, it
would be desirable to find simplified analytical expressions when the special
functions are evaluated on arguments with values close to unity. Second, as
demonstrated by the example of the triple-exchange channel provided in
Appendix D and for which we derived preliminary results, our technique can be
used to obtain fully analytical results for diagrams involving more-than-two
exchanges of massive fields, definitely providing fascinating directions for
future work. Third, it would be intriguing to investigate whether the well-
known Suyama–Yamaguchi inequality [190], relating the collapsed limit of the
trispectrum to the squeezed one of the bispectrum, still holds in the context
of heavy field exchanges, depending on the number of such exchanges. Fourth,
as already acknowledged in the main text and following the now standard
tradition in computing cosmological correlators, we focused here on the
exchange of massive fields in pure dS, but future works on multiple exchange
diagrams should also include symmetry-breaking cases, such as those with non-
unit sound speed, chemical potential, IR effects, and time-dependent masses
and coupling constants. Fifth, although we did mention inflationary flavor
oscillations in the context of double-exchange channels and uncovered first
interesting unique features, it remains yet to conduct a thorough
investigation of how to tell different channels apart in practice from (to-be)
observations. Sixth, with the analytical results of tree-level double-exchange
here, and the known result of the 1-loop bubble diagram [136] that also
contains two massive propagators, it would be interesting to compare the
signals originating from these different channels. Seventh, we explained in
the phenomenology Section 5 that the cubic interaction $\lambda$, relevant for
double-exchange channels, can be related to the field-space curvature in
general non-linear sigma models of inflation; it would be thrilling to explore
double-exchange correlators in the context of concrete models of moduli field
spaces in string compactifications or the coset field space of (pseudo-)Nambu-
Goldstones from spontaneous symmetry breaking patterns. Lastly, only very
recently was the CC signal arising from varied single-exchange diagrams
searched for in the cosmological data [156, 157], but given the relatively
higher prospects of detectability of the double-exchange channels as we
explained at length, it seems important to extend the search to templates
proposed in this work. We plan to address several of these directions in
future works.

##### Acknowledgements.

We would like to thank Kohei Fujikura, Toshifumi Noumi, Zhehan Qin, Sébastien
Renaux-Petel, David Stefanyszyn, Xi Tong, Denis Werth and Zhong-Zhi Xianyu for
useful discussion. L.P. would like to acknowledge Denis Werth for helpful
discussions regarding the CosmoFlow code used in this work. S.A. is supported
by IBS under the project code, IBS-R018-D1. L.P. acknowledges funding support
from the Initiative Physique des Infinis (IPI), a research training program of
the Idex SUPER at Sorbonne Université. F.S. acknowledges financial aid from
the Institute for Basic Science, and JSPS Grant-in-Aid for Scientific Research
No. 23KJ0938. M.Y. is supported by IBS under the project code, IBS-R018-D3,
and by JSPS Grant-in-Aid for Scientific Research Number JP21H01080. Y.Z. is
supported by the IBS under the project code, IBS-R018-D3.

##  Appendix A Mathematical Formulae

In this Appendix, we provide a summary of definitions and formulae frequently
utilized in this work. Most of these formulae can also be found in
Mathematical functions handbook [173].

Following Refs. [135, 138, 141], we use shorthand notations for products of
Gamma function,

$\displaystyle\Gamma\left[z_{1},\cdots,z_{m}\right]\equiv\Gamma\left(z_{1}\right)\cdots\Gamma\left(z_{m}\right),$
(A.1) $\displaystyle\Gamma\left[\begin{array}[]{c}z_{1},\cdots,z_{m}\\\
w_{1},\cdots,w_{n}\end{array}\right]\equiv\frac{\Gamma\left(z_{1}\right)\cdots\Gamma\left(z_{m}\right)}{\Gamma\left(w_{1}\right)\cdots\Gamma\left(w_{n}\right)}.$
(A.4)

Several types of hypergeometric functions are used in the main text, which we
collect their definitions here. First, the (generalized) hypergeometric
function ${}_{p}\rm{F}_{q}$ is defined by

$\displaystyle{}_{p}{\rm{F}}_{q}\left[\begin{array}[]{c|c}a_{1},\cdots,a_{p}\\\
b_{1},\cdots,b_{q}\end{array}\
x\right]=\sum_{n=0}^{\infty}\frac{\left(a_{1}\right)_{n}\cdots\left(a_{p}\right)_{n}}{\left(b_{1}\right)_{n}\cdots\left(b_{q}\right)_{n}}\frac{x^{n}}{n!},$
(A.7)

where $(a)_{n}\equiv\Gamma[a+n]/\Gamma[a]$ is the Pochhammer symbol.

We also use its “dressed” version, which is defined by

$\displaystyle{}_{p}\mathcal{F}_{q}\left[\begin{array}[]{c|c}a_{1},\cdots,a_{p}\\\
b_{1},\cdots,b_{q}\end{array}\ x\right]$
$\displaystyle\equiv\Gamma\left[\begin{array}[]{c}a_{1},\cdots,a_{p}\\\
b_{1},\cdots,b_{q}\end{array}\right]{}_{p}\mathrm{F}_{q}\left[\begin{array}[]{c|c}a_{1},\cdots,a_{p}\\\
b_{1},\cdots,b_{q}\end{array}\ x\right]$ (A.14)
$\displaystyle=\sum_{n=0}^{\infty}\Gamma\left[\begin{array}[]{c}a_{1}+n,\cdots,a_{p}+n\\\
b_{1}+n,\cdots,b_{q}+n\end{array}\right]\frac{x^{n}}{n!}.$ (A.17)

There are many formulae changing the variables of hypergeometric functions.
The one frequently used in the main text is

$\displaystyle\frac{\sin{(\pi(c-a-b))}}{\pi}{}_{2}\mathcal{F}_{1}\left[\begin{array}[]{c|c}a,b\\\
c\end{array}\ x\right]$ (A.20)
$\displaystyle=\frac{1}{\Gamma\left[c-a,c-b\right]}\,{}_{2}\mathcal{F}_{1}\left[\begin{array}[]{c|c}a,b\\\
a+b-c+1\end{array}\
1-x\right]-\frac{(1-x)^{c-a-b}}{\Gamma\left[c-a,c-b\right]}{}_{2}\mathcal{F}_{1}\left[\begin{array}[]{c|c}c-a,c-b\\\
c-a-b+1\end{array}\ 1-x\right]~{}.$ (A.25)

Next, let us move to the hypergeometric series of two variables. In this
paper, we used two types of Appell series, $F_{2}$ and $F_{4}$, whose
definitions are given by

$\displaystyle F_{2}\left[a\left|\begin{array}[]{l}b_{1},b_{2}\\\
c_{1},c_{2}\end{array}\right|x,y\right]$
$\displaystyle=\sum_{m,n=0}^{\infty}\frac{(a)_{m+n}\left(b_{1}\right)_{m}\left(b_{2}\right)_{n}}{\left(c_{1}\right)_{m}\left(c_{2}\right)_{n}m!n!}x^{m}y^{n},$
(A.28) $\displaystyle F_{4}\left[\begin{array}[]{c|c}a,b\\\
c_{1},c_{2}\end{array}\ x,y\right]$
$\displaystyle=\sum_{m,n=0}^{\infty}\frac{(a)_{m+n}(b)_{m+n}}{\left(c_{1}\right)_{m}\left(c_{2}\right)_{n}m!n!}x^{m}y^{n},$
(A.31)

and their dressed versions are expressed as

$\displaystyle\mathcal{F}_{2}\left[a\left|\begin{array}[]{l}b_{1},b_{2}\\\
c_{1},c_{2}\end{array}\right|x,y\right]$
$\displaystyle=\sum_{m,n=0}^{\infty}\Gamma\left[\begin{array}[]{c}a+m+n,b_{1}+m,b_{2}+n\\\
c_{1}+m,c_{2}+n\end{array}\right]\frac{x^{m}y^{n}}{m!n!},$ (A.36)
$\displaystyle\mathcal{F}_{4}\left[\begin{array}[]{c|c}a,b\\\
c_{1},c_{2}\end{array}\ x,y\right]$
$\displaystyle=\sum_{m,n=0}^{\infty}\Gamma\left[\begin{array}[]{c}a+m+n,b+m+n\\\
c_{1}+m,c_{2}+n\end{array}\right]\frac{x^{m}y^{n}}{m!n!}.$ (A.41)

Note that the function $F_{2}$ is convergent for $|x|+|y|<1$, while the
function $F_{4}$ for $\sqrt{|x|}+\sqrt{|y|}<1$. Here are some useful
conversion formulae. The Appell $F_{4}$ function is related to Appell $F_{2}$
through

$\displaystyle F_{4}\left[\begin{array}[]{c|c}\frac{a}{2},\frac{a+1}{2}\\\
b_{1}+\frac{1}{2},b_{2}+\frac{1}{2}\end{array}\
x^{2},y^{2}\right]=(1+x+y)^{-a}F_{2}\left[a\left|\begin{array}[]{l}b_{1},b_{2}\\\
2b_{1},2b_{2}\end{array}\right|\frac{2x}{x+y+1},\frac{2y}{x+y+1}\right]~{},$
(A.46)

and the Appell $F_{2}$ can be expanded as summation of Gauss hypergeometric
function ${}_{2}{\rm{F}}_{1}$ as

$\displaystyle F_{2}\left[a\left|\begin{array}[]{l}b_{1},b_{2}\\\
c_{1},c_{2}\end{array}\right|x,y\right]=\sum_{m=0}^{\infty}\frac{(a)_{m}(b_{1})_{m}}{(c_{1})_{m}}\frac{x^{m}}{m!}{}_{2}{\rm{F}}_{1}\left[\begin{array}[]{c|c}a+m,b_{2}\\\
c_{2}\end{array}\ y\right]~{}.$ (A.51)

Finally, the $F_{C}$-type of Lauricella’s hypergeometric functions of $m$
variables $x_{1},\cdots,x_{m}$ appears in Appendix D, which is defined by
[191],

$\displaystyle F_{C}\left[\begin{array}[]{c|c}a,b\\\
c_{1},\cdots,c_{m}\end{array}\
x_{1},\cdots,x_{m}\right]=\sum_{n_{1},\ldots,n_{m}=0}^{\infty}\frac{(a)_{n_{1}+\cdots+n_{m}}(b)_{n_{1}+\cdots+n_{m}}}{\left(c_{1}\right)_{n_{1}}\cdots\left(c_{m}\right)_{n_{m}}n_{1}!\cdots
n_{m}!}x_{1}^{n_{1}}\cdots x_{m}^{n_{m}},$ (A.54)

with convergent radius

$\displaystyle\sqrt{\left|x_{1}\right|}+\cdots+\sqrt{\left|x_{m}\right|}<1.$
(A.55)

Note that Lauricella’s $F_{C}$ function is reduced to Appell’s and Gauss
hypergeometric functions $F_{4}$ and ${}_{2}{\rm{F}}_{1}$ for $m=2$ and $m=1$,
respectively.

##  Appendix B Seed Integral with Single Massive Exchange

Here we show analytic expressions of the seed integral with single exchange,
which appear as source terms of bootstrap equations for double-exchange: the
right-hand side of (2.51), (2.53), and (2.54). The analytical solution of the
single-exchange tree-level diagram is well-established in the literature, as
we mentioned in the introduction. Further, in cases where one vertex involves
linear mixing, the solution provides an exact closed form. This solution as
the source terms is utilized extensively in our main text discussion. The
comprehensive derivation can be found in detail in [138]. Here, we just
collect the results in a convenient format for our purpose.

To do so, let us transform the seed integral with single exchange (2.44) by
changing the variables

$\displaystyle z_{1}\equiv-k_{124}\tau_{1},\quad z_{2}=-k_{3}\tau_{2},$ (B.1)

and introducing

$\displaystyle
D_{\mathsf{ab}}^{\alpha}\left(k_{3};\tau_{1},\tau_{2}\right)=k_{3}^{-3}\widehat{D}^{\alpha}_{\mathsf{ab}}\left(\frac{R}{2-R}z_{1},z_{2}\right),$
(B.2)

where

$\displaystyle R\equiv\frac{2k_{3}}{k_{1234}}.$ (B.3)

Under these change, Eq. (2.44) can be written as

$\displaystyle\mathcal{I}_{\mathsf{ab},\alpha}^{p_{1}p_{2}}(R)$
$\displaystyle=H^{-2}\left(\frac{R}{2-R}\right)^{1+p_{1}}(-\mathsf{ab})\int_{0}^{\infty}\mathrm{d}z_{1}\mathrm{d}z_{2}\
z_{1}^{p_{1}}z_{2}^{p_{2}}e^{-\mathrm{i}\mathsf{a}z_{1}-\mathrm{i}\mathsf{b}z_{2}}\widehat{D}^{\alpha}_{\mathsf{ab}}\left(\frac{R}{2-R}z_{1},z_{2}\right).$
(B.4)

Note that the seed integral depends only on a specific momentum ratio $R$.

A summarized version of the explicit expression for Eq. (B.4) as a function of
$R$ is given by

$\displaystyle\mathcal{I}_{\pm\mp,\alpha}^{p_{1}p_{2}}(R)=$ $\displaystyle\
\frac{e^{\mp\mathrm{i}\frac{\pi}{2}\bar{p}_{12}}}{2^{6+p_{12}-2\mathrm{i}\mu_{\alpha}}\pi^{\frac{1}{2}}}\Gamma\left[\begin{array}[]{c}\frac{5}{2}+p_{2}+\mathrm{i}\mu_{\alpha},\frac{5}{2}+p_{2}-\mathrm{i}\mu_{\alpha},\frac{5}{2}+p_{1}-\mathrm{i}\mu_{\alpha},\mathrm{i}\mu_{\alpha}\\\
3+p_{2}\end{array}\right]$ (B.7) $\displaystyle\times
R^{\frac{5}{2}+p_{1}-\mathrm{i}\mu_{\alpha}}{}_{2}{\rm{F}}_{1}\left[\begin{array}[]{c|c}\frac{1}{2}-\mathrm{i}\mu_{\alpha},\frac{5}{2}+p_{1}-\mathrm{i}\mu_{\alpha}\\\
1-2\mathrm{i}\mu_{\alpha}\end{array}\
R\right]+(\mu_{\alpha}\rightarrow-\mu_{\alpha}),$ (B.10)
$\displaystyle\mathcal{I}_{\pm\pm,\alpha}^{p_{1}p_{2}}(R)=$
$\displaystyle\left\\{\frac{\pm\mathrm{i}\
e^{\mp\mathrm{i}\frac{\pi}{2}p_{12}}e^{\mp\pi\mu_{\alpha}}}{2^{6+p_{12}-2\mathrm{i}\mu_{\alpha}}\pi^{\frac{1}{2}}}\Gamma\left[\begin{array}[]{c}\frac{5}{2}+p_{2}+\mathrm{i}\mu_{\alpha},\frac{5}{2}+p_{2}-\mathrm{i}\mu_{\alpha},\frac{5}{2}+p_{1}-\mathrm{i}\mu_{\alpha},\mathrm{i}\mu_{\alpha}\\\
3+p_{2}\end{array}\right]\right.$ (B.13) $\displaystyle\left.\times
R^{\frac{5}{2}+p_{1}-\mathrm{i}\mu_{\alpha}}{}_{2}{\rm{F}}_{1}\left[\begin{array}[]{c|c}\frac{1}{2}-\mathrm{i}\mu_{\alpha},\frac{5}{2}+p_{1}-\mathrm{i}\mu_{\alpha}\\\
1-2\mathrm{i}\mu_{\alpha}\end{array}\
R\right]+(\mu_{\alpha}\rightarrow-\mu_{\alpha})\right\\}$ (B.16)
$\displaystyle+\frac{e^{\mp\mathrm{i}\frac{\pi}{2}p_{12}}\Gamma\left(p_{12}+5\right)}{\mu_{\alpha}^{2}+\left(\frac{5}{2}+p_{2}\right)^{2}}\left(\frac{R}{2}\right)^{5+p_{12}}{}_{3}{\rm{F}}_{2}\left[\begin{array}[]{c|c}1,5+p_{12},3+p_{2}\\\
\frac{7}{2}+p_{2}-\mathrm{i}\mu_{\alpha},\frac{7}{2}+p_{2}+\mathrm{i}\mu_{\alpha}\end{array}\
R\right],$ (B.19)

where $p_{12}\equiv p_{1}+p_{2}$ and $\bar{p}_{12}\equiv p_{1}-p_{2}$.

##  Appendix C Boundary Conditions of Seed Integrals

In the main text, in Section 3.2, we imposed the soft limit ($u,v\rightarrow
0$) as the boundary conditions to determine the coefficients of bootstrap
equations. In this Appendix, we present the calculation of the time integrals
under these limits. In principle, one can employ the late-time expansion of
the Hankel function to directly evaluate the time integral under these limits.
Instead, as suggested by previous work [135], the partial Mellin-Barnes (MB)
representation proves to be an efficient approach for handling such time
integrals. So in this Appendix we will adopt the MB representation to find the
boundary condition. We will not extensively discuss the derivation details,
but instead, we will present the final results in a format suitable for our
purposes. Further details on the use of the partial MB method to evaluate
double-exchange correlators can be found in [141] .
The key trick of this method is to use the following representation of Hankel
function,

$\displaystyle\mathrm{H}_{\mathrm{i}\mu}^{(j)}(-k\tau)=\int_{-\mathrm{i}\infty}^{\mathrm{i}\infty}\frac{\mathrm{d}s}{2\pi\mathrm{i}}\frac{(-k\tau/2)^{-2s}}{\pi}e^{(-1)^{j+1}(2s-\mathrm{i}\mu-1)\mathrm{i}\frac{\pi}{2}}\Gamma\left[s-\frac{\mathrm{i}\mu}{2},s+\frac{\mathrm{i}\mu}{2}\right],\quad(j=1,2)$
(C.1)

by which the time-integral in the seed integral can be trivially performed.
Then, the remaining complex integral with variable $s$ can be done by the
residue theorem. With (C.1), the factorised propagators in the seed integral
(2.10) can be expressed as

$\displaystyle D^{\alpha}_{\pm\mp}\left(k_{1};\tau_{1},\tau_{2}\right)=$
$\displaystyle\ \frac{H^{2}\pi
e^{-\pi\mu_{\alpha}}}{4}\left(\tau_{1}\tau_{2}\right)^{3/2}H_{\mp\mathrm{i}\mu_{\alpha}}^{(2,1)}\left(-k_{1}\tau_{1}\right)H_{\pm\mathrm{i}\mu_{\alpha}}^{(1,2)}\left(-k_{1}\tau_{2}\right)$
$\displaystyle=$ $\displaystyle\
\frac{H^{2}}{4\pi}\int_{-\mathrm{i}\infty}^{\mathrm{i}\infty}\frac{\mathrm{d}s_{1}}{2\pi\mathrm{i}}\frac{\mathrm{d}s_{2}}{2\pi\mathrm{i}}e^{\mp\mathrm{i}\pi\left(s_{1}-s_{2}\right)}\left(\frac{k_{1}}{2}\right)^{-2s_{12}}\left(-\tau_{1}\right)^{-2s_{1}+3/2}\left(-\tau_{2}\right)^{-2s_{2}+3/2}$
$\displaystyle\times\Gamma\left[s_{1}-\frac{\mathrm{i}\mu_{\alpha}}{2},s_{1}+\frac{\mathrm{i}\mu_{\alpha}}{2},s_{2}-\frac{\mathrm{i}\mu_{\alpha}}{2},s_{2}+\frac{\mathrm{i}\mu_{\alpha}}{2}\right],$
(C.2)

and similarly for $D^{\beta}_{\pm\mp}\left(k_{3};\tau_{2},\tau_{3}\right)$.
For time-ordered propagators, it has been observed that a more effective way
involves dividing them into factorised and time-ordered components through

$\displaystyle
D^{\alpha}_{\pm\pm}\left(k_{1};\tau_{1},\tau_{2}\right)=D^{\alpha}_{\pm\mp}\left(k_{1};\tau_{1},\tau_{2}\right)+\left[D^{\alpha}_{\mp\pm}\left(k_{1};\tau_{1},\tau_{2}\right)-D^{\alpha}_{\pm\mp}\left(k_{1};\tau_{1},\tau_{2}\right)\right]\theta\left(\tau_{1}-\tau_{2}\right),$
(C.3) $\displaystyle
D^{\beta}_{\pm\pm}\left(k_{3};\tau_{2},\tau_{3}\right)=D^{\beta}_{\mp\pm}\left(k_{3};\tau_{2},\tau_{3}\right)+\left[D^{\beta}_{\pm\mp}\left(k_{3};\tau_{2},\tau_{3}\right)-D^{\beta}_{\mp\pm}\left(k;\tau_{2},\tau_{3}\right)\right]\theta\left(\tau_{3}-\tau_{2}\right).$
(C.4)

Then accordingly, we define factorised and time-ordered parts of the seed
integral by [141]

$\displaystyle\mathcal{I}_{\pm\mp\pm}~{},$ (C.5)
$\displaystyle\mathcal{I}_{\pm\pm\mp}=\mathcal{I}_{\pm\pm\mp,\mathrm{F}}+\mathcal{I}_{\pm\pm\mp,\mathrm{T}}~{},$
(C.6)
$\displaystyle\mathcal{I}_{\pm\mp\mp}=\mathcal{I}_{\pm\mp\mp,\mathrm{F}}+\mathcal{I}_{\pm\mp\mp,\mathrm{T}}~{},$
(C.7)
$\displaystyle\mathcal{I}_{\pm\pm\pm}=\mathcal{I}_{\pm\pm\pm,\mathrm{FF}}+\mathcal{I}_{\pm\pm\pm,\mathrm{TF}}+\mathcal{I}_{\pm\pm\pm,\mathrm{FT}}+\mathcal{I}_{\pm\pm\pm,\mathrm{TT}}~{},$
(C.8)

where $\mathcal{I}_{\pm\mp\pm}$ is already factorised from the beginning. The
subscript F(T) indicates that it originates from the factorised (time-ordered)
part of the time-ordered propagators. In the soft limit, the (fully)
factorised part will dominate [141]. To clarify further, the dominant parts
are

$\displaystyle\mathcal{I}_{\pm\mp\pm}=$
$\displaystyle\pm\mathrm{i}H^{-4}k_{24}^{9+p_{123}}\int^{0}_{-\infty}\mathrm{d}\tau_{1}\mathrm{d}\tau_{2}\mathrm{d}\tau_{3}(-\tau_{1})^{p_{1}}(-\tau_{2})^{p_{2}}(-\tau_{3})^{p_{3}}e^{\pm\mathrm{i}k_{1}\tau_{1}\mp\mathrm{i}k_{24}\tau_{2}\pm\mathrm{i}k_{3}\tau_{3}}$
$\displaystyle\times
D^{\alpha}_{\pm\mp}\left(k_{1};\tau_{1},\tau_{2}\right)D^{\beta}_{\mp\pm}\left(k_{3};\tau_{2},\tau_{3}\right),$
(C.9) $\displaystyle\mathcal{I}_{\pm\pm\mp,\mathrm{F}}=$
$\displaystyle\pm\mathrm{i}H^{-4}k_{24}^{9+p_{123}}\int^{0}_{-\infty}\mathrm{d}\tau_{1}\mathrm{d}\tau_{2}\mathrm{d}\tau_{3}(-\tau_{1})^{p_{1}}(-\tau_{2})^{p_{2}}(-\tau_{3})^{p_{3}}e^{\pm\mathrm{i}k_{1}\tau_{1}\pm\mathrm{i}k_{24}\tau_{2}\mp\mathrm{i}k_{3}\tau_{3}}$
$\displaystyle\times
D^{\alpha}_{\pm\mp}\left(k_{1};\tau_{1},\tau_{2}\right)D^{\beta}_{\pm\mp}\left(k_{3};\tau_{2},\tau_{3}\right),$
(C.10) $\displaystyle\mathcal{I}_{\pm\mp\mp,\mathrm{F}}=$
$\displaystyle\mp\mathrm{i}H^{-4}k_{24}^{9+p_{123}}\int^{0}_{-\infty}\mathrm{d}\tau_{1}\mathrm{d}\tau_{2}\mathrm{d}\tau_{3}(-\tau_{1})^{p_{1}}(-\tau_{2})^{p_{2}}(-\tau_{3})^{p_{3}}e^{\pm\mathrm{i}k_{1}\tau_{1}\mp\mathrm{i}k_{24}\tau_{2}\mp\mathrm{i}k_{3}\tau_{3}}$
$\displaystyle\times
D^{\alpha}_{\pm\mp}\left(k_{1};\tau_{1},\tau_{2}\right)D^{\beta}_{\pm\mp}\left(k_{3};\tau_{2},\tau_{3}\right),$
(C.11) $\displaystyle\mathcal{I}_{\pm\pm\pm,\mathrm{FF}}=$
$\displaystyle\mp\mathrm{i}H^{-4}k_{24}^{9+p_{123}}\int^{0}_{-\infty}\mathrm{d}\tau_{1}\mathrm{d}\tau_{2}\mathrm{d}\tau_{3}(-\tau_{1})^{p_{1}}(-\tau_{2})^{p_{2}}(-\tau_{3})^{p_{3}}e^{\pm\mathrm{i}k_{1}\tau_{1}\pm\mathrm{i}k_{24}\tau_{2}\pm\mathrm{i}k_{3}\tau_{3}}$
$\displaystyle\times
D^{\alpha}_{\pm\mp}\left(k_{1};\tau_{1},\tau_{2}\right)D^{\beta}_{\mp\pm}\left(k_{3};\tau_{2},\tau_{3}\right).$
(C.12)

As mentioned above, the things to do are to insert the MB representations of
the SK propagators (C.2) into the seed integrals above, perform $\tau_{i}\
(i=1,2,3)$-integral, and evaluate the complex $s_{j}\ (j=1,2,3,4)$-integral by
residue theorem. After straightforward calculation, we gets the results
summarized below.

As for the fully factorised seed integral
$\mathcal{I}_{\pm\mp\pm,\alpha\beta}^{p_{1}p_{2}p_{3}}$, the final results are

$\displaystyle\mathcal{I}_{\pm\mp\pm,\alpha\beta}^{p_{1}p_{2}p_{3}}=$
$\displaystyle\
\frac{e^{\mp\mathrm{i}\frac{\pi}{2}(p_{13}-p_{2})}}{2^{7+p_{13}}\pi}u^{-p_{1}-\frac{5}{2}}v^{-p_{3}-\frac{5}{2}}\Gamma\left[\begin{array}[]{c}\frac{5}{2}+p_{1}-\mathrm{i}\mu_{\alpha},\frac{5}{2}+p_{1}+\mathrm{i}\mu_{\alpha},\frac{5}{2}+p_{3}-\mathrm{i}\mu_{\beta},\frac{5}{2}+p_{3}+\mathrm{i}\mu_{\beta}\\\
3+p_{1},3+p_{3}\end{array}\right]$ (C.15)
$\displaystyle\times\sum_{\mathsf{a},\mathsf{b}=\pm}\
\sum_{n_{2},n_{3}=0}^{\infty}\left(\frac{u}{2}\right)^{2n_{2}-\mathrm{i}\mathsf{a}\mu_{\alpha}}\left(\frac{v}{2}\right)^{2n_{3}-\mathrm{i}\mathsf{b}\mu_{\beta}}\frac{(-1)^{n_{23}}}{n_{2}!n_{3}!}$
$\displaystyle\times\Gamma\left[-n_{2}+\mathrm{i}\mathsf{a}\mu_{\alpha},-n_{3}+\mathrm{i}\mathsf{b}\mu_{\beta},p_{2}+4+2n_{23}-\mathrm{i}(\mathsf{a}\mu_{\alpha}+\mathsf{b}\mu_{\beta})\right],$
(C.16)

where $n_{23}=n_{2}+n_{3}$. Under the limit $u,v\ll 1$ for boundary conditions
of bootstrap equations, the dominant contribution comes from $n_{2,3}=0$ in
summations. Then, utilizing the following formulae for the Gamma function,

$\displaystyle\Gamma[z,1-z]=\frac{\pi}{\sin\pi
z},\quad\Gamma\left[\begin{array}[]{c}2z\\\
z,z+\frac{1}{2}\end{array}\right]=\frac{2^{2z-1}}{\pi^{1/2}},$ (C.19)

one may summarize

$\displaystyle\lim_{u,v\ll
1}\mathcal{I}_{\pm\mp\pm,\alpha\beta}^{p_{1}p_{2}p_{3}}=$
$\displaystyle-\sum_{\mathsf{a,b}=\pm}e^{\mp\mathrm{i}\frac{\pi}{2}\left(p_{13}-p_{2}\right)}\operatorname{csch}\left(\pi\mathsf{a}\mu_{\alpha}\right)\operatorname{csch}\left(\pi\mathsf{b}\mu_{\beta}\right)\widetilde{\Gamma}(p_{1},p_{2},p_{3},\mu_{\alpha},\mu_{\beta})$
$\displaystyle\times
u^{-\frac{5}{2}-p_{1}-\mathrm{i}\mathsf{a}\mu_{\alpha}}v^{-\frac{5}{2}-p_{3}-\mathrm{i}\mathsf{b}\mu_{\beta}}\
\Gamma\left[\begin{array}[]{c}\frac{4+p_{2}-\mathrm{i}\left(\mathsf{a}\mu_{\alpha}+\mathsf{b}\mu_{\beta}\right)}{2},\frac{5+p_{2}-\mathrm{i}\left(\mathsf{a}\mu_{\alpha}+\mathsf{b}\mu_{\beta}\right)}{2}\\\
1-\mathrm{i}\mathsf{a}\mu_{\alpha},1-\mathrm{i}\mathsf{b}\mu_{\beta}\end{array}\right],$
(C.22)

with

$\displaystyle\widetilde{\Gamma}(p_{1},p_{2},p_{3},\mu_{\alpha},\mu_{\beta})\equiv\frac{\pi^{\frac{1}{2}}}{2^{4+p_{13}-p_{2}}}\Gamma\left[\begin{array}[]{c}\frac{5}{2}+p_{1}-\mathrm{i}\mu_{\alpha},\frac{5}{2}+p_{1}+\mathrm{i}\mu_{\alpha},\frac{5}{2}+p_{3}-\mathrm{i}\mu_{\beta},\frac{5}{2}+p_{3}+\mathrm{i}\mu_{\beta}\\\
3+p_{1},3+p_{3}\end{array}\right]~{}.$ (C.25)

A similar procedure can be applied to both the partially factorised (partially
nested) and fully nested seed integrals. Under double soft limits, they yield:

$\displaystyle\lim_{u,v\ll
1}\mathcal{I}_{\pm\pm\mp,\alpha\beta,}^{p_{1}p_{2}p_{3}}=$
$\displaystyle\mp\mathrm{i}\sum_{\mathsf{a,b}=\pm}e^{\mp\mathrm{i}\frac{\pi}{2}\left(p_{12}-p_{3}\right)}\operatorname{csch}\left(\pi\mathsf{a}\mu_{\alpha}\right)\operatorname{csch}\left(\pi\mathsf{b}\mu_{\beta}\right)e^{\mp\pi\mathsf{a}\mu_{\alpha}}\widetilde{\Gamma}(p_{1},p_{2},p_{3},\mu_{\alpha},\mu_{\beta})$
$\displaystyle\times
u^{-\frac{5}{2}-p_{1}-\mathrm{i}\mathsf{a}\mu_{\alpha}}v^{-\frac{5}{2}-p_{3}-\mathrm{i}\mathsf{b}\mu_{\beta}}\
\Gamma\left[\begin{array}[]{c}\frac{4+p_{2}-\mathrm{i}\left(\mathsf{a}\mu_{\alpha}+\mathsf{b}\mu_{\beta}\right)}{2},\frac{5+p_{2}-\mathrm{i}\left(\mathsf{a}\mu_{\alpha}+\mathsf{b}\mu_{\beta}\right)}{2}\\\
1-\mathrm{i}\mathsf{a}\mu_{\alpha},1-\mathrm{i}\mathsf{b}\mu_{\beta}\end{array}\right],$
(C.28) $\displaystyle\lim_{u,v\ll
1}\mathcal{I}_{\pm\pm\pm,\alpha\beta}^{p_{1}p_{2}p_{3}}=$
$\displaystyle\sum_{\mathsf{a,b}=\pm}e^{\mp\mathrm{i}\frac{\pi}{2}p_{123}}\operatorname{csch}\left(\pi\mathsf{a}\mu_{\alpha}\right)\operatorname{csch}\left(\pi\mathsf{b}\mu_{\beta}\right)e^{\mp\pi(\mathsf{a}\mu_{\alpha}+\mathsf{b}\mu_{\beta})}\widetilde{\Gamma}(p_{1},p_{2},p_{3},\mu_{\alpha},\mu_{\beta})$
$\displaystyle\times
u^{-\frac{5}{2}-p_{1}-\mathrm{i}\mathsf{a}\mu_{\alpha}}v^{-\frac{5}{2}-p_{3}-\mathrm{i}\mathsf{b}\mu_{\beta}}\
\Gamma\left[\begin{array}[]{c}\frac{4+p_{2}-\mathrm{i}\left(\mathsf{a}\mu_{\alpha}+\mathsf{b}\mu_{\beta}\right)}{2},\frac{5+p_{2}-\mathrm{i}\left(\mathsf{a}\mu_{\alpha}+\mathsf{b}\mu_{\beta}\right)}{2}\\\
1-\mathrm{i}\mathsf{a}\mu_{\alpha},1-\mathrm{i}\mathsf{b}\mu_{\beta}\end{array}\right].$
(C.31)

The results for $\mathcal{I}_{\pm\mp\mp,\alpha\beta}^{p_{1}p_{2}p_{3}}$ can be
directly obtained by replacing $u\leftrightarrow v$, $p_{1}\leftrightarrow
p_{3}$, and $\alpha\leftrightarrow\beta$ in the expression of
$\mathcal{I}_{\mp\mp\pm,\alpha\beta}^{p_{1}p_{2}p_{3}}$.

##  Appendix D Toward Triple Exchange

In principle, the bootstrap equation approach we applied to the double-
exchange diagram in this work can be easily generalized to any tree-level
diagrams. In this section, let us see how it works for the triple-exchange
diagram, keeping in mind that obtaining the complete analytic form remains a
non-trivial task for future work. In particular, the bispectrum with triple
massive exchange would be phenomenologically important because it can give a
sizable signal known, as first emphasized in the context of quasi-single-field
inflation [1].

Here again, we start from the seed integral for the four-point inflaton
correlator, but from a triple exchange, which can be defined by

$\displaystyle\mathcal{I}_{\mathsf{abcd},\alpha\beta\gamma}^{p_{1}p_{2}p_{3}p_{4}}=$
$\displaystyle\ H^{-6}k_{4}^{13+p_{1234}}(\mathsf{abcd})\int^{0}_{-\infty}\
\prod_{i=1}^{4}\ \mathrm{d}\tau_{i}\
(-\tau_{i})^{p_{i}}e^{\mathrm{i}\mathsf{a}k_{1}\tau_{1}+\mathrm{i}\mathsf{b}k_{2}\tau_{2}+\mathrm{i}\mathsf{c}k_{3}\tau_{3}+\mathrm{i}\mathsf{d}k_{4}\tau_{4}}$
$\displaystyle\times
D^{\alpha}_{\mathsf{ad}}\left(k_{1};\tau_{1},\tau_{4}\right)D^{\beta}_{\mathsf{bd}}\left(k_{2};\tau_{2},\tau_{4}\right)D^{\gamma}_{\mathsf{cd}}\left(k_{3};\tau_{3},\tau_{4}\right).$
(D.1)

It can be diagrammatically expressed as shown in Figure 12. The three-point
correlator with triple exchange can be obtained by setting the external
momentum soft, i.e., $k_{4}\rightarrow 0$. In the following, we consider a
situation where three massive fields contribute to the seed from the mixing
vertex $\sigma^{\alpha}\sigma^{\beta}\sigma^{\gamma}$. Without ambiguity and
for the sake of brevity, we will omit mass indices ($\alpha,\beta,\gamma$) in
the following text.

Figure 12: Four-point and three-point diagrams with triple massive exchange.
For generality, distinct colors here denote massive scalar fields
$\sigma^{\alpha}$, $\sigma^{\beta}$ and $\sigma^{\gamma}$ with different mass.

Similarly to the procedure shown in the main text, we transform the seed (D.1)
to a more convenient form through changing the integration variables by

$\displaystyle-k_{i}\tau_{i}=z_{i},\quad(i=1,2,3,4),$ (D.2)

and define the hat propagator in the same way as Eqs. (2.22)-(2.24). Then, the
seed integral (D.1) can be written as

$\displaystyle\mathcal{I}_{\mathsf{abcd}}^{p_{1}p_{2}p_{3}p_{4}}=$
$\displaystyle\
\frac{H^{-6}(\mathsf{abcd})}{w_{1}^{4+p_{1}}w_{2}^{4+p_{2}}w_{3}^{4+p_{3}}}\int_{0}^{\infty}\prod_{i=1}^{4}\
\mathrm{d}z_{i}\
z_{i}^{p_{i}}e^{-\mathrm{i}\mathsf{a}z_{1}-\mathrm{i}\mathsf{b}z_{2}-\mathrm{i}\mathsf{c}z_{3}-\mathrm{i}\mathsf{d}z_{4}}\widehat{D}_{\mathsf{ad}}\left(z_{1},w_{1}z_{4}\right)\widehat{D}_{\mathsf{bd}}\left(z_{2},w_{2}z_{4}\right)\widehat{D}_{\mathsf{cd}}\left(z_{3},w_{3}z_{4}\right),$
$\displaystyle\equiv$ $\displaystyle\
\frac{1}{w_{1}^{4+p_{1}}w_{2}^{4+p_{2}}w_{3}^{4+p_{3}}}\widehat{\mathcal{I}}_{\mathsf{abcd}}^{p_{1}p_{2}p_{3}p_{4}}(w_{1},w_{2},w_{3})~{},$
(D.3)

where

$\displaystyle w_{i}\equiv\frac{k_{i}}{k_{4}}~{},\quad(i=1,2,3).$ (D.4)

Now let us derive the bootstrap equations for $\widehat{\mathcal{I}}$. For
example, repeating the same procedure as in the main text, we find
$\widehat{\mathcal{I}}_{\pm\pm\pm\mp}$ to satisfy the following homogeneous
differential equations

$\displaystyle\mathcal{D}^{(3)}_{w_{i}}\,\widehat{\mathcal{I}}_{\pm\pm\pm\mp}^{p_{1}p_{2}p_{3}p_{4}}=0,~{}~{}\qquad(i=1,2,3)~{},$
(D.5)

where the differential operators $\mathcal{D}^{(3)}_{w_{i}}\ (i=1,2,3)$ are
defined by

$\displaystyle\mathcal{D}^{(3)}_{w_{i}}$ $\displaystyle\equiv
w_{i}^{2}\partial_{w_{i}}^{2}-2w_{i}\partial_{w_{i}}+\mu_{i}^{2}+\frac{9}{4}-w_{i}^{2}\left(\sum_{j=1}^{3}w_{j}\partial_{w_{j}}+p_{4}+2\right)\left(\sum_{j=1}^{3}w_{j}\partial_{w_{j}}+p_{4}+1\right)$
$\displaystyle=w_{i}^{2}\left(1-w_{i}^{2}\right)\partial_{w_{i}}^{2}-2w_{i}^{3}\partial_{w_{i}}\sum_{j\neq
i}w_{j}\partial_{w_{j}}-w_{i}^{2}\sum_{j,k\neq
i}w_{j}w_{k}\partial_{w_{j}}\partial_{w_{k}}-2w_{i}\left(1+\left(p_{4}+2\right)w_{i}^{2}\right)\partial_{w_{i}}$
$\displaystyle\ \ \ \ -2\left(p_{4}+2\right)w_{i}^{2}\sum_{j\neq
i}w_{j}\partial_{w_{j}}+\mu_{i}^{2}+\frac{9}{4}-(p_{4}+2)(p_{4}+1)w_{i}^{2}\,,$
(D.6)

where $\mu_{i}$ is the mass parameter of the massive field associated with the
momentum $k_{i}$. Remarkably, this can be analytically solved by

$\displaystyle\widehat{\mathcal{I}}_{\pm\pm\pm\mp}^{p_{1}p_{2}p_{3}p_{4}}=$
$\displaystyle\ \sum_{\mathsf{abc}=\pm}c_{\pm\pm\pm\mp,\mathsf{abc}}\
w_{1}^{\frac{3}{2}-\mathrm{i}\mathsf{a}\mu_{1}}w_{2}^{\frac{3}{2}-\mathrm{i}\mathsf{b}\mu_{2}}w_{3}^{\frac{3}{2}-\mathrm{i}\mathsf{c}\mu_{3}}$
$\displaystyle\times
F_{C}\left[\begin{array}[]{c|c}\frac{2p_{4}+11-2\mathrm{i}(\mathsf{a}\mu_{1}+\mathsf{b}\mu_{2}+\mathsf{c}\mu_{3})}{4},\frac{2p_{4}+13-2\mathrm{i}(\mathsf{a}\mu_{1}+\mathsf{b}\mu_{2}+\mathsf{c}\mu_{3})}{4}\\\
1-\mathrm{i}\mathsf{a}\mu_{1},1-\mathrm{i}\mathsf{b}\mu_{2},1-\mathrm{i}\mathsf{c}\mu_{3}\end{array}\
w_{1}^{2},w_{2}^{2},w_{3}^{2}\right],$ (D.9)

where the function $F_{C}[\cdots]$ is Lauricella’s hypergeometric function
with three variables defined in (A.54). The eight coefficients
$c_{\pm\pm\pm\mp,\mathsf{abc}}$ could be determined by proper boundary
conditions.

Figure 13: The schematic diagram illustrates the triple-exchange bootstrap
equations. $\mathsf{a},\mathsf{b},\mathsf{c}$ and $\mathsf{d}$ are the SK
indices. When the two SK indices of one massive propagator are opposite,
applying the corresponding differential operator results in homogeneous
equations. Conversely, if the two indices are the same, an additional
$\delta$-term introduces a source term, which can be expressed using the seed
integral of the double-exchange diagram.

As a next example, let us see $\widehat{\mathcal{I}}_{\mp\pm\pm\mp}$ including
a single time-ordered propagator for $\sigma^{\alpha}$,
$\widehat{D}_{\mp\mp}(z_{1},r_{1}z_{4})$. In this case, we obtain an
inhomogeneous term in the bootstrap equation after acting with the
differential operator $\mathcal{D}^{(3)}_{w_{1}}$,

$\displaystyle\mathcal{D}^{(3)}_{w_{1}}\,\widehat{\mathcal{I}}_{\mp\pm\pm\mp,\alpha\beta\gamma}^{p_{1}p_{2}p_{3}p_{4}}=\frac{w_{1}^{4+p_{1}}w_{2}^{4+p_{2}}w_{3}^{4+p_{3}}}{\left(1+w_{1}\right)^{p_{1234}+13}}\mathcal{I}_{\pm\mp\pm}^{p_{2},p_{14}+4,p_{3}}\left(u=\frac{w_{2}}{1+w_{1}},v=\frac{w_{3}}{1+w_{1}}\right)~{},$
(D.10)

while the other equations with $\mathcal{D}^{(3)}_{w_{2}}$ and
$\mathcal{D}^{(3)}_{w_{3}}$ still satisfy the homogeneous ones,

$\displaystyle\mathcal{D}^{(3)}_{w_{2}}\,\widehat{\mathcal{I}}_{\mp\pm\pm\mp,\alpha\beta\gamma}^{p_{1}p_{2}p_{3}p_{4}}=\mathcal{D}^{(3)}_{w_{3}}\,\widehat{\mathcal{I}}_{\mp\pm\pm\mp,\alpha\beta\gamma}^{p_{1}p_{2}p_{3}p_{4}}=0~{}.$
(D.11)

Note that the right-hand side of the first equation (D.10) is written as the
seed integral of double-exchange diagram, which we obtained analytic forms in
the main text. This structure is generic: the seed integral with $n$-massive
field exchange, denoted as $\mathcal{I}_{n}$, emerges as a source term in the
bootstrap equations of $\mathcal{I}_{n+1}$. Consequently, we can
systematically solve more complex tree-level diagrams, although the detailed
calculations remain non-trivial tasks: as shown in (D.9), the homogeneous
solutions can be written down as known special functions, which is true even
for more general diagram with $n$-massive exchanges. However, determining
particular solutions would be a non-trivial task since the source term (the
right-hand side of (D.10) for triple exchange) becomes more and more
complicated as $n$ increases, and we leave the detailed derivation for future
work.
Similarly, deriving the remaining bootstrap equations for other seed integrals
follows a straightforward procedure. However, we do not intend to delve into
all the tedious details within this work. Instead, we illustrate the structure
of the bootstrap equations in Figure 13.

## References

  * [1] X. Chen and Y. Wang, _Quasi-Single Field Inflation and Non-Gaussianities_ , _JCAP_ 04 (2010) 027 [0911.3380].
  * [2] D. Baumann and D. Green, _Signatures of Supersymmetry from the Early Universe_ , _Phys. Rev. D_ 85 (2012) 103520 [1109.0292].
  * [3] T. Noumi, M. Yamaguchi and D. Yokoyama, _Effective field theory approach to quasi-single field inflation and effects of heavy fields_ , _JHEP_ 06 (2013) 051 [1211.1624].
  * [4] N. Arkani-Hamed and J. Maldacena, _Cosmological Collider Physics_ , 1503.08043.
  * [5] X. Chen and Y. Wang, _Large non-Gaussianities with Intermediate Shapes from Quasi-Single Field Inflation_ , _Phys. Rev. D_ 81 (2010) 063511 [0909.0496].
  * [6] E. Sefusatti, J.R. Fergusson, X. Chen and E.P.S. Shellard, _Effects and Detectability of Quasi-Single Field Inflation in the Large-Scale Structure and Cosmic Microwave Background_ , _JCAP_ 08 (2012) 033 [1204.6318].
  * [7] J. Norena, L. Verde, G. Barenboim and C. Bosch, _Prospects for constraining the shape of non-Gaussianity with the scale-dependent bias_ , _JCAP_ 08 (2012) 019 [1204.6324].
  * [8] X. Chen and Y. Wang, _Quasi-Single Field Inflation with Large Mass_ , _JCAP_ 09 (2012) 021 [1205.0160].
  * [9] S. Pi and M. Sasaki, _Curvature Perturbation Spectrum in Two-field Inflation with a Turning Trajectory_ , _JCAP_ 10 (2012) 051 [1205.0161].
  * [10] S. Céspedes and G.A. Palma, _Cosmic inflation in a landscape of heavy-fields_ , _JCAP_ 10 (2013) 051 [1303.4703].
  * [11] J.-O. Gong, S. Pi and M. Sasaki, _Equilateral non-Gaussianity from heavy fields_ , _JCAP_ 11 (2013) 043 [1306.3691].
  * [12] R. Emami, _Spectroscopy of Masses and Couplings during Inflation_ , _JCAP_ 04 (2014) 031 [1311.0184].
  * [13] A. Kehagias and A. Riotto, _High Energy Physics Signatures from Inflation and Conformal Symmetry of de Sitter_ , _Fortsch. Phys._ 63 (2015) 531 [1501.03515].
  * [14] J. Liu, Y. Wang and S. Zhou, _Inflation with Massive Vector Fields_ , _JCAP_ 08 (2015) 033 [1502.05138].
  * [15] E. Dimastrogiovanni, M. Fasiello and M. Kamionkowski, _Imprints of Massive Primordial Fields on Large-Scale Structure_ , _JCAP_ 02 (2016) 017 [1504.05993].
  * [16] F. Schmidt, N.E. Chisari and C. Dvorkin, _Imprint of inflation on galaxy shape correlations_ , _JCAP_ 10 (2015) 032 [1506.02671].
  * [17] X. Chen, M.H. Namjoo and Y. Wang, _Quantum Primordial Standard Clocks_ , _JCAP_ 02 (2016) 013 [1509.03930].
  * [18] L.V. Delacretaz, T. Noumi and L. Senatore, _Boost Breaking in the EFT of Inflation_ , _JCAP_ 02 (2017) 034 [1512.04100].
  * [19] B. Bonga, S. Brahma, A.-S. Deutsch and S. Shandera, _Cosmic variance in inflation with two light scalars_ , _JCAP_ 05 (2016) 018 [1512.05365].
  * [20] X. Chen, Y. Wang and Z.-Z. Xianyu, _Loop Corrections to Standard Model Fields in Inflation_ , _JHEP_ 08 (2016) 051 [1604.07841].
  * [21] R. Flauger, M. Mirbabayi, L. Senatore and E. Silverstein, _Productive Interactions: heavy particles and non-Gaussianity_ , _JCAP_ 10 (2017) 058 [1606.00513].
  * [22] H. Lee, D. Baumann and G.L. Pimentel, _Non-Gaussianity as a Particle Detector_ , _JHEP_ 12 (2016) 040 [1607.03735].
  * [23] L.V. Delacretaz, V. Gorbenko and L. Senatore, _The Supersymmetric Effective Field Theory of Inflation_ , _JHEP_ 03 (2017) 063 [1610.04227].
  * [24] P.D. Meerburg, M. Münchmeyer, J.B. Muñoz and X. Chen, _Prospects for Cosmological Collider Physics_ , _JCAP_ 03 (2017) 050 [1610.06559].
  * [25] X. Chen, Y. Wang and Z.-Z. Xianyu, _Standard Model Background of the Cosmological Collider_ , _Phys. Rev. Lett._ 118 (2017) 261302 [1610.06597].
  * [26] X. Chen, Y. Wang and Z.-Z. Xianyu, _Standard Model Mass Spectrum in Inflationary Universe_ , _JHEP_ 04 (2017) 058 [1612.08122].
  * [27] A. Kehagias and A. Riotto, _On the Inflationary Perturbations of Massive Higher-Spin Fields_ , _JCAP_ 07 (2017) 046 [1705.05834].
  * [28] H. An, M. McAneny, A.K. Ridgway and M.B. Wise, _Quasi Single Field Inflation in the non-perturbative regime_ , _JHEP_ 06 (2018) 105 [1706.09971].
  * [29] X. Tong, Y. Wang and S. Zhou, _On the Effective Field Theory for Quasi-Single Field Inflation_ , _JCAP_ 11 (2017) 045 [1708.01709].
  * [30] A. Moradinezhad Dizgah and C. Dvorkin, _Scale-Dependent Galaxy Bias from Massive Particles with Spin during Inflation_ , _JCAP_ 01 (2018) 010 [1708.06473].
  * [31] A.V. Iyer, S. Pi, Y. Wang, Z. Wang and S. Zhou, _Strongly Coupled Quasi-Single Field Inflation_ , _JCAP_ 01 (2018) 041 [1710.03054].
  * [32] H. An, M. McAneny, A.K. Ridgway and M.B. Wise, _Non-Gaussian Enhancements of Galactic Halo Correlations in Quasi-Single Field Inflation_ , _Phys. Rev. D_ 97 (2018) 123528 [1711.02667].
  * [33] S. Kumar and R. Sundrum, _Heavy-Lifting of Gauge Theories By Cosmic Inflation_ , _JHEP_ 05 (2018) 011 [1711.03988].
  * [34] G. Franciolini, A. Kehagias and A. Riotto, _Imprints of Spinning Particles on Primordial Cosmological Perturbations_ , _JCAP_ 02 (2018) 023 [1712.06626].
  * [35] S. Riquelme M., _Non-Gaussianities in a two-field generalization of Natural Inflation_ , _JCAP_ 04 (2018) 027 [1711.08549].
  * [36] X. Tong, Y. Wang and S. Zhou, _Unsuppressed primordial standard clocks in warm quasi-single field inflation_ , _JCAP_ 06 (2018) 013 [1801.05688].
  * [37] X. Chen, W.Z. Chua, Y. Guo, Y. Wang, Z.-Z. Xianyu and T. Xie, _Quantum Standard Clocks in the Primordial Trispectrum_ , _JCAP_ 05 (2018) 049 [1803.04412].
  * [38] R. Saito and T. Kubota, _Heavy Particle Signatures in Cosmological Correlation Functions with Tensor Modes_ , _JCAP_ 06 (2018) 009 [1804.06974].
  * [39] G. Cabass, E. Pajer and F. Schmidt, _Imprints of Oscillatory Bispectra on Galaxy Clustering_ , _JCAP_ 09 (2018) 003 [1804.07295].
  * [40] A. Moradinezhad Dizgah, H. Lee, J.B. Muñoz and C. Dvorkin, _Galaxy Bispectrum from Massive Spinning Particles_ , _JCAP_ 05 (2018) 013 [1801.07265].
  * [41] Y. Wang, Y.-P. Wu, J. Yokoyama and S. Zhou, _Hybrid Quasi-Single Field Inflation_ , _JCAP_ 07 (2018) 068 [1804.07541].
  * [42] X. Chen, Y. Wang and Z.-Z. Xianyu, _Neutrino Signatures in Primordial Non-Gaussianities_ , _JHEP_ 09 (2018) 022 [1805.02656].
  * [43] N. Bartolo, D.M. Bianco, R. Jimenez, S. Matarrese and L. Verde, _Supergravity, $\alpha$-attractors and primordial non-Gaussianity_, _JCAP_ 10 (2018) 017 [1805.04269].
  * [44] E. Dimastrogiovanni, M. Fasiello and G. Tasinato, _Probing the inflationary particle content: extra spin-2 field_ , _JCAP_ 08 (2018) 016 [1806.00850].
  * [45] L. Bordin, P. Creminelli, A. Khmelnitsky and L. Senatore, _Light Particles with Spin in Inflation_ , _JCAP_ 10 (2018) 013 [1806.10587].
  * [46] X. Chen, A. Loeb and Z.-Z. Xianyu, _Unique Fingerprints of Alternatives to Inflation in the Primordial Power Spectrum_ , _Phys. Rev. Lett._ 122 (2019) 121301 [1809.02603].
  * [47] A. Achúcarro, S. Céspedes, A.-C. Davis and G.A. Palma, _Constraints on Holographic Multifield Inflation and Models Based on the Hamilton-Jacobi Formalism_ , _Phys. Rev. Lett._ 122 (2019) 191301 [1809.05341].
  * [48] W.Z. Chua, Q. Ding, Y. Wang and S. Zhou, _Imprints of Schwinger Effect on Primordial Spectra_ , _JHEP_ 04 (2019) 066 [1810.09815].
  * [49] S. Kumar and R. Sundrum, _Seeing Higher-Dimensional Grand Unification In Primordial Non-Gaussianities_ , _JHEP_ 04 (2019) 120 [1811.11200].
  * [50] G. Goon, K. Hinterbichler, A. Joyce and M. Trodden, _Shapes of gravity: Tensor non-Gaussianity and massive spin-2 fields_ , _JHEP_ 10 (2019) 182 [1812.07571].
  * [51] Y.-P. Wu, _Higgs as heavy-lifted physics during inflation_ , _JHEP_ 04 (2019) 125 [1812.10654].
  * [52] D. Anninos, V. De Luca, G. Franciolini, A. Kehagias and A. Riotto, _Cosmological Shapes of Higher-Spin Gravity_ , _JCAP_ 04 (2019) 045 [1902.01251].
  * [53] L. Li, T. Nakama, C.M. Sou, Y. Wang and S. Zhou, _Gravitational Production of Superheavy Dark Matter and Associated Cosmological Signatures_ , _JHEP_ 07 (2019) 067 [1903.08842].
  * [54] M. McAneny and A.K. Ridgway, _New Shapes of Primordial Non-Gaussianity from Quasi-Single Field Inflation with Multiple Isocurvatons_ , _Phys. Rev. D_ 100 (2019) 043534 [1903.11607].
  * [55] S. Kim, T. Noumi, K. Takeuchi and S. Zhou, _Heavy Spinning Particles from Signs of Primordial Non-Gaussianities: Beyond the Positivity Bounds_ , _JHEP_ 12 (2019) 107 [1906.11840].
  * [56] S. Alexander, S.J. Gates, L. Jenks, K. Koutrolikos and E. McDonough, _Higher Spin Supersymmetry at the Cosmological Collider: Sculpting SUSY Rilles in the CMB_ , _JHEP_ 10 (2019) 156 [1907.05829].
  * [57] S. Lu, Y. Wang and Z.-Z. Xianyu, _A Cosmological Higgs Collider_ , _JHEP_ 02 (2020) 011 [1907.07390].
  * [58] A. Hook, J. Huang and D. Racco, _Searches for other vacua. Part II. A new Higgstory at the cosmological collider_ , _JHEP_ 01 (2020) 105 [1907.10624].
  * [59] A. Hook, J. Huang and D. Racco, _Minimal signatures of the Standard Model in non-Gaussianities_ , _Phys. Rev. D_ 101 (2020) 023519 [1908.00019].
  * [60] S. Kumar and R. Sundrum, _Cosmological Collider Physics and the Curvaton_ , _JHEP_ 04 (2020) 077 [1908.11378].
  * [61] T. Liu, X. Tong, Y. Wang and Z.-Z. Xianyu, _Probing P and CP Violations on the Cosmological Collider_ , _JHEP_ 04 (2020) 189 [1909.01819].
  * [62] L.-T. Wang and Z.-Z. Xianyu, _In Search of Large Signals at the Cosmological Collider_ , _JHEP_ 02 (2020) 044 [1910.12876].
  * [63] Y. Welling, _Simple, exact model of quasisingle field inflation_ , _Phys. Rev. D_ 101 (2020) 063535 [1907.02951].
  * [64] D.-G. Wang, _On the inflationary massive field with a curved field manifold_ , _JCAP_ 01 (2020) 046 [1911.04459].
  * [65] Y. Wang and Y. Zhu, _Cosmological Collider Signatures of Massive Vectors from Non-Gaussian Gravitational Waves_ , _JCAP_ 04 (2020) 049 [2001.03879].
  * [66] L. Li, S. Lu, Y. Wang and S. Zhou, _Cosmological Signatures of Superheavy Dark Matter_ , _JHEP_ 07 (2020) 231 [2002.01131].
  * [67] L.-T. Wang and Z.-Z. Xianyu, _Gauge Boson Signals at the Cosmological Collider_ , _JHEP_ 11 (2020) 082 [2004.02887].
  * [68] J. Fan and Z.-Z. Xianyu, _A Cosmic Microscope for the Preheating Era_ , _JHEP_ 01 (2021) 021 [2005.12278].
  * [69] K. Kogai, K. Akitsu, F. Schmidt and Y. Urakawa, _Galaxy imaging surveys as spin-sensitive detector for cosmological colliders_ , _JCAP_ 03 (2021) 060 [2009.05517].
  * [70] A. Bodas, S. Kumar and R. Sundrum, _The Scalar Chemical Potential in Cosmological Collider Physics_ , _JHEP_ 02 (2021) 079 [2010.04727].
  * [71] S. Aoki and M. Yamaguchi, _Disentangling mass spectra of multiple fields in cosmological collider_ , _JHEP_ 04 (2021) 127 [2012.13667].
  * [72] N. Maru and A. Okawa, _Non-Gaussianity from $X,Y$ gauge bosons in Cosmological Collider Physics_, 2101.10634.
  * [73] S. Kim, T. Noumi, K. Takeuchi and S. Zhou, _Perturbative unitarity in quasi-single field inflation_ , _JHEP_ 07 (2021) 018 [2102.04101].
  * [74] S. Lu, _Axion isocurvature collider_ , _JHEP_ 04 (2022) 157 [2103.05958].
  * [75] C.M. Sou, X. Tong and Y. Wang, _Chemical-potential-assisted particle production in FRW spacetimes_ , _JHEP_ 06 (2021) 129 [2104.08772].
  * [76] Q. Lu, M. Reece and Z.-Z. Xianyu, _Missing scalars at the cosmological collider_ , _JHEP_ 12 (2021) 098 [2108.11385].
  * [77] L. Pinol, S. Aoki, S. Renaux-Petel and M. Yamaguchi, _Inflationary flavor oscillations and the cosmic spectroscopy_ , _Phys. Rev. D_ 107 (2023) L021301 [2112.05710].
  * [78] Y. Cui and Z.-Z. Xianyu, _Probing Leptogenesis with the Cosmological Collider_ , _Phys. Rev. Lett._ 129 (2022) 111301 [2112.10793].
  * [79] X. Tong and Z.-Z. Xianyu, _Large spin-2 signals at the cosmological collider_ , _JHEP_ 10 (2022) 194 [2203.06349].
  * [80] M. Reece, L.-T. Wang and Z.-Z. Xianyu, _Large-Field Inflation and the Cosmological Collider_ , 2204.11869.
  * [81] X. Chen, R. Ebadi and S. Kumar, _Classical cosmological collider physics and primordial features_ , _JCAP_ 08 (2022) 083 [2205.01107].
  * [82] Z. Qin and Z.-Z. Xianyu, _Phase information in cosmological collider signals_ , _JHEP_ 10 (2022) 192 [2205.01692].
  * [83] G. Cabass, S. Jazayeri, E. Pajer and D. Stefanyszyn, _Parity violation in the scalar trispectrum: no-go theorems and yes-go examples_ , 2210.02907.
  * [84] G. Cabass, M.M. Ivanov and O.H.E. Philcox, _Colliding Ghosts: Constraining Inflation with the Parity-Odd Galaxy Four-Point Function_ , 2210.16320.
  * [85] X. Niu, M.H. Rahat, K. Srinivasan and W. Xue, _Gravitational Wave Probes of Massive Gauge Bosons at the Cosmological Collider_ , 2211.14331.
  * [86] X. Niu, M.H. Rahat, K. Srinivasan and W. Xue, _Parity-Odd and Even Trispectrum from Axion Inflation_ , 2211.14324.
  * [87] D.H. Tran, Y. Wang, J. Yang and Y. Zhu, _Numerical prescriptions of early-time divergences of the in-in formalism_ , _JCAP_ 07 (2022) 014 [2202.06350].
  * [88] S. Aoki, _Continuous spectrum on cosmological collider_ , _JCAP_ 04 (2023) 002 [2301.07920].
  * [89] X. Tong, Y. Wang, C. Zhang and Y. Zhu, _BCS in the Sky: Signatures of Inflationary Fermion Condensation_ , 2304.09428.
  * [90] S. Jazayeri, S. Renaux-Petel and D. Werth, _Shapes of the cosmological low-speed collider_ , _JCAP_ 12 (2023) 035 [2307.01751].
  * [91] Y. Yin, _The cosmological collider signal in the non-BD initial states_ , 2309.05244.
  * [92] P. Chakraborty and J. Stout, _Light Scalars at the Cosmological Collider_ , 2310.01494.
  * [93] X. Chen, J. Fan and L. Li, _New inflationary probes of axion dark matter_ , 2303.03406.
  * [94] S. Jazayeri, S. Renaux-Petel, X. Tong, D. Werth and Y. Zhu, _Parity violation from emergent nonlocality during inflation_ , _Phys. Rev. D_ 108 (2023) 123523 [2308.11315].
  * [95] P. Chakraborty and J. Stout, _Compact Scalars at the Cosmological Collider_ , 2311.09219.
  * [96] Y. Ema and S. Verner, _Cosmological Collider Signatures of Higgs- $R^{2}$ Inflation_, 2309.10841.
  * [97] C. McCulloch, E. Pajer and X. Tong, _A Cosmological Tachyon Collider: Enhancing the Long-Short Scale Coupling_ , 2401.11009.
  * [98] N. Craig, S. Kumar and A. McCune, _An Effective Cosmological Collider_ , 2401.10976.
  * [99] Y.-P. Wu, _The cosmological collider in $R^{2}$ inflation_, 2404.05031.
  * [100] J.M. Maldacena and G.L. Pimentel, _On graviton non-Gaussianities during inflation_ , _JHEP_ 09 (2011) 045 [1104.2846].
  * [101] V. Assassi, D. Baumann and D. Green, _On Soft Limits of Inflationary Correlation Functions_ , _JCAP_ 11 (2012) 047 [1204.4207].
  * [102] D. Baumann, G. Goon, H. Lee and G.L. Pimentel, _Partially Massless Fields During Inflation_ , _JHEP_ 04 (2018) 140 [1712.06624].
  * [103] N. Arkani-Hamed and P. Benincasa, _On the Emergence of Lorentz Invariance and Unitarity from the Scattering Facet of Cosmological Polytopes_ , 1811.01125.
  * [104] C. Sleight, _A Mellin Space Approach to Cosmological Correlators_ , _JHEP_ 01 (2020) 090 [1906.12302].
  * [105] C. Sleight and M. Taronna, _Bootstrapping Inflationary Correlators in Mellin Space_ , _JHEP_ 02 (2020) 098 [1907.01143].
  * [106] N. Arkani-Hamed, D. Baumann, H. Lee and G.L. Pimentel, _The Cosmological Bootstrap: Inflationary Correlators from Symmetries and Singularities_ , _JHEP_ 04 (2020) 105 [1811.00024].
  * [107] D. Baumann, C. Duaso Pueyo, A. Joyce, H. Lee and G.L. Pimentel, _The cosmological bootstrap: weight-shifting operators and scalar seeds_ , _JHEP_ 12 (2020) 204 [1910.14051].
  * [108] D. Baumann, C. Duaso Pueyo, A. Joyce, H. Lee and G.L. Pimentel, _The Cosmological Bootstrap: Spinning Correlators from Symmetries and Factorization_ , _SciPost Phys._ 11 (2021) 071 [2005.04234].
  * [109] E. Pajer, D. Stefanyszyn and J. Supeł, _The Boostless Bootstrap: Amplitudes without Lorentz boosts_ , _JHEP_ 12 (2020) 198 [2007.00027].
  * [110] H. Goodhew, S. Jazayeri and E. Pajer, _The Cosmological Optical Theorem_ , _JCAP_ 04 (2021) 021 [2009.02898].
  * [111] D. Meltzer, _The inflationary wavefunction from analyticity and factorization_ , _JCAP_ 12 (2021) 018 [2107.10266].
  * [112] H. Goodhew, S. Jazayeri, M.H.G. Lee and E. Pajer, _Cutting cosmological correlators_ , _JCAP_ 08 (2021) 003 [2104.06587].
  * [113] S. Jazayeri, E. Pajer and D. Stefanyszyn, _From locality and unitarity to cosmological correlators_ , _JHEP_ 10 (2021) 065 [2103.08649].
  * [114] D. Baumann, W.-M. Chen, C. Duaso Pueyo, A. Joyce, H. Lee and G.L. Pimentel, _Linking the singularities of cosmological correlators_ , _JHEP_ 09 (2022) 010 [2106.05294].
  * [115] S. Melville and E. Pajer, _Cosmological Cutting Rules_ , _JHEP_ 05 (2021) 249 [2103.09832].
  * [116] L. Di Pietro, V. Gorbenko and S. Komatsu, _Analyticity and unitarity for cosmological correlators_ , _JHEP_ 03 (2022) 023 [2108.01695].
  * [117] H. Gomez, R.L. Jusinskas and A. Lipstein, _Cosmological Scattering Equations_ , _Phys. Rev. Lett._ 127 (2021) 251604 [2106.11903].
  * [118] A. Premkumar, _Regulating loops in de Sitter spacetime_ , _Phys. Rev. D_ 109 (2024) 045003 [2110.12504].
  * [119] J. Bonifacio, E. Pajer and D.-G. Wang, _From amplitudes to contact cosmological correlators_ , _JHEP_ 10 (2021) 001 [2106.15468].
  * [120] M. Hogervorst, J.a. Penedones and K.S. Vaziri, _Towards the non-perturbative cosmological bootstrap_ , 2107.13871.
  * [121] C. Sleight and M. Taronna, _From AdS to dS exchanges: Spectral representation, Mellin amplitudes, and crossing_ , _Phys. Rev. D_ 104 (2021) L081902 [2007.09993].
  * [122] C. Sleight and M. Taronna, _From dS to AdS and back_ , _JHEP_ 12 (2021) 074 [2109.02725].
  * [123] C. Sleight and M. Taronna, _On the consistency of (partially-)massless matter couplings in de Sitter space_ , _JHEP_ 10 (2021) 156 [2106.00366].
  * [124] G. Cabass, E. Pajer, D. Stefanyszyn and J. Supeł, _Bootstrapping large graviton non-Gaussianities_ , _JHEP_ 05 (2022) 077 [2109.10189].
  * [125] X. Tong, Y. Wang and Y. Zhu, _Cutting rule for cosmological collider signals: a bulk evolution perspective_ , _JHEP_ 03 (2022) 181 [2112.03448].
  * [126] L.-T. Wang, Z.-Z. Xianyu and Y.-M. Zhong, _Precision calculation of inflation correlators at one loop_ , _JHEP_ 02 (2022) 085 [2109.14635].
  * [127] H. Gomez, R. Lipinski Jusinskas and A. Lipstein, _Cosmological scattering equations at tree-level and one-loop_ , _JHEP_ 07 (2022) 004 [2112.12695].
  * [128] D. Baumann, D. Green, A. Joyce, E. Pajer, G.L. Pimentel, C. Sleight et al., _Snowmass White Paper: The Cosmological Bootstrap_ , in _Snowmass 2021_ , 3, 2022 [2203.08121].
  * [129] T. Heckelbacher, I. Sachs, E. Skvortsov and P. Vanhove, _Analytical evaluation of cosmological correlation functions_ , _JHEP_ 08 (2022) 139 [2204.07217].
  * [130] E. Pajer, _Building a Boostless Bootstrap for the Bispectrum_ , _JCAP_ 01 (2021) 023 [2010.12818].
  * [131] G.L. Pimentel and D.-G. Wang, _Boostless cosmological collider bootstrap_ , _JHEP_ 10 (2022) 177 [2205.00013].
  * [132] S. Jazayeri and S. Renaux-Petel, _Cosmological bootstrap in slow motion_ , _JHEP_ 12 (2022) 137 [2205.10340].
  * [133] J. Bonifacio, H. Goodhew, A. Joyce, E. Pajer and D. Stefanyszyn, _The graviton four-point function in de Sitter space_ , _JHEP_ 06 (2023) 212 [2212.07370].
  * [134] H. Lee and X. Wang, _Cosmological double-copy relations_ , _Phys. Rev. D_ 108 (2023) L061702 [2212.11282].
  * [135] Z. Qin and Z.-Z. Xianyu, _Helical inflation correlators: partial Mellin-Barnes and bootstrap equations_ , _JHEP_ 04 (2023) 059 [2208.13790].
  * [136] Z.-Z. Xianyu and H. Zhang, _Bootstrapping one-loop inflation correlators with the spectral decomposition_ , _JHEP_ 04 (2023) 103 [2211.03810].
  * [137] D.-G. Wang, G.L. Pimentel and A. Achúcarro, _Bootstrapping multi-field inflation: non-Gaussianities from light scalars revisited_ , _JCAP_ 05 (2023) 043 [2212.14035].
  * [138] Z. Qin and Z.-Z. Xianyu, _Closed-form formulae for inflation correlators_ , _JHEP_ 07 (2023) 001 [2301.07047].
  * [139] Z. Qin and Z.-Z. Xianyu, _Inflation correlators at the one-loop order: nonanalyticity, factorization, cutting rule, and OPE_ , _JHEP_ 09 (2023) 116 [2304.13295].
  * [140] Z. Qin and Z.-Z. Xianyu, _Nonanalyticity and on-shell factorization of inflation correlators at all loop orders_ , _JHEP_ 01 (2024) 168 [2308.14802].
  * [141] Z.-Z. Xianyu and J. Zang, _Inflation Correlators with Multiple Massive Exchanges_ , 2309.10849.
  * [142] D. Green, Y. Huang, C.-H. Shen and D. Baumann, _Positivity from Cosmological Correlators_ , 2310.02490.
  * [143] C. Duaso Pueyo and E. Pajer, _A Cosmological Bootstrap for Resonant Non-Gaussianity_ , 2311.01395.
  * [144] S. De and A. Pokraka, _Cosmology meets cohomology_ , 2308.03753.
  * [145] M. Loparco, J. Penedones, K. Salehi Vaziri and Z. Sun, _The Källén-Lehmann representation in de Sitter spacetime_ , _JHEP_ 12 (2023) 159 [2306.00090].
  * [146] M.H.G. Lee, C. McCulloch and E. Pajer, _Leading loops in cosmological correlators_ , _JHEP_ 11 (2023) 038 [2305.11228].
  * [147] S. Aoki, T. Noumi, F. Sano and M. Yamaguchi, _Analytic Formulae for Inflationary Correlators with Dynamical Mass_ , 2312.09642.
  * [148] N. Arkani-Hamed, D. Baumann, A. Hillman, A. Joyce, H. Lee and G.L. Pimentel, _Differential Equations for Cosmological Correlators_ , 2312.05303.
  * [149] N. Arkani-Hamed, D. Baumann, A. Hillman, A. Joyce, H. Lee and G.L. Pimentel, _Kinematic Flow and the Emergence of Time_ , 2312.05300.
  * [150] D. Stefanyszyn, X. Tong and Y. Zhu, _Cosmological Correlators Through the Looking Glass: Reality, Parity, and Factorisation_ , 2309.07769.
  * [151] Y. Donath and E. Pajer, _The In-Out Formalism for In-In Correlators_ , 2402.05999.
  * [152] B. Fan and Z.-Z. Xianyu, _Cosmological Amplitudes in Power-Law FRW Universe_ , 2403.07050.
  * [153] T.W. Grimm, A. Hoefnagels and M. van Vliet, _Structure and Complexity of Cosmological Correlators_ , 2404.03716.
  * [154] Z. Du and D. Stefanyszyn, _Soft Theorems for Boostless Amplitudes_ , 2403.05459.
  * [155] S. Melville and G.L. Pimentel, _A de Sitter S-matrix from amputated cosmological correlators_ , 2404.05712.
  * [156] G. Cabass, O.H.E. Philcox, M.M. Ivanov, K. Akitsu, S.-F. Chen, M. Simonović et al., _BOSS Constraints on Massive Particles during Inflation: The Cosmological Collider in Action_ , 2404.01894.
  * [157] W. Sohn, D.-G. Wang, J.R. Fergusson and E.P.S. Shellard, _Searching for Cosmological Collider in the Planck CMB Data_ , 2404.07203.
  * [158] X. Chen, _Primordial Non-Gaussianities from Inflation Models_ , _Adv. Astron._ 2010 (2010) 638979 [1002.1416].
  * [159] Y. Wang, _Inflation, Cosmic Perturbations and Non-Gaussianities_ , _Commun. Theor. Phys._ 62 (2014) 109 [1303.1523].
  * [160] X. Chen, Y. Wang and Z.-Z. Xianyu, _Schwinger-Keldysh Diagrammatics for Primordial Perturbations_ , _JCAP_ 12 (2017) 006 [1703.10166].
  * [161] D. Anninos, T. Anous, D.Z. Freedman and G. Konstantinidis, _Late-time Structure of the Bunch-Davies De Sitter Wavefunction_ , _JCAP_ 11 (2015) 048 [1406.5490].
  * [162] N. Arkani-Hamed, P. Benincasa and A. Postnikov, _Cosmological Polytopes and the Wavefunction of the Universe_ , 1709.02813.
  * [163] A. Hillman, _Symbol Recursion for the dS Wave Function_ , 1912.09450.
  * [164] S. Céspedes, A.-C. Davis and S. Melville, _On the time evolution of cosmological correlators_ , _JHEP_ 02 (2021) 012 [2009.07874].
  * [165] A. Hillman and E. Pajer, _A differential representation of cosmological wavefunctions_ , _JHEP_ 04 (2022) 012 [2112.01619].
  * [166] S.A. Salcedo, M.H.G. Lee, S. Melville and E. Pajer, _The Analytic Wavefunction_ , _JHEP_ 06 (2023) 020 [2212.08009].
  * [167] S. Albayrak, P. Benincasa and C. Duaso Pueyo, _Perturbative Unitarity and the Wavefunction of the Universe_ , 2305.19686.
  * [168] S. Agui-Salcedo and S. Melville, _The cosmological tree theorem_ , _JHEP_ 12 (2023) 076 [2308.00680].
  * [169] S. Céspedes, A.-C. Davis and D.-G. Wang, _On the IR Divergences in de Sitter Space: loops, resummation and the semi-classical wavefunction_ , 2311.17990.
  * [170] D. Werth, L. Pinol and S. Renaux-Petel, _Cosmological Flow of Primordial Correlators_ , 2302.00655.
  * [171] L. Pinol, S. Renaux-Petel and D. Werth, _The Cosmological Flow: A Systematic Approach to Primordial Correlators_ , 2312.06559.
  * [172] D. Werth, L. Pinol and S. Renaux-Petel, _CosmoFlow: Python Package for Cosmological Correlators_ , 2402.03693.
  * [173] “NIST Digital Library of Mathematical Functions.” https://dlmf.nist.gov/.
  * [174] S. Garcia-Saenz, L. Pinol and S. Renaux-Petel, _Revisiting non-Gaussianity in multifield inflation with curved field space_ , _JHEP_ 01 (2020) 073 [1907.10403].
  * [175] L. Pinol, _Multifield inflation beyond $N_{\mathrm{field}}=2$: non-Gaussianities and single-field effective theory_, _JCAP_ 04 (2021) 002 [2011.05930].
  * [176] A. Bzowski, P. McFadden and K. Skenderis, _Implications of conformal invariance in momentum space_ , _JHEP_ 03 (2014) 111 [1304.7760].
  * [177] W. Bühring, _Generalized hypergeometric functions at unit argument_ , _Proceedings of the American Mathematical Society_ 114 (1992) 145.
  * [178] R. Holman and A.J. Tolley, _Enhanced Non-Gaussianity from Excited Initial States_ , _JCAP_ 05 (2008) 001 [0710.1302].
  * [179] R. Flauger, D. Green and R.A. Porto, _On squeezed limits in single-field inflation. Part I_ , _JCAP_ 08 (2013) 032 [1303.1430].
  * [180] Planck collaboration, _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10 [1807.06211].
  * [181] P. Creminelli, M.A. Luty, A. Nicolis and L. Senatore, _Starting the Universe: Stable Violation of the Null Energy Condition and Non-standard Cosmologies_ , _JHEP_ 12 (2006) 080 [hep-th/0606090].
  * [182] C. Cheung, P. Creminelli, A.L. Fitzpatrick, J. Kaplan and L. Senatore, _The Effective Field Theory of Inflation_ , _JHEP_ 03 (2008) 014 [0709.0293].
  * [183] L. Senatore and M. Zaldarriaga, _The Effective Field Theory of Multifield Inflation_ , _JHEP_ 04 (2012) 024 [1009.2093].
  * [184] A.J. Tolley and M. Wyman, _The Gelaton Scenario: Equilateral non-Gaussianity from multi-field dynamics_ , _Phys. Rev. D_ 81 (2010) 043502 [0910.1853].
  * [185] S. Cremonini, Z. Lalak and K. Turzynski, _Strongly Coupled Perturbations in Two-Field Inflationary Models_ , _JCAP_ 03 (2011) 016 [1010.3021].
  * [186] A. Achucarro, J.-O. Gong, S. Hardeman, G.A. Palma and S.P. Patil, _Features of heavy physics in the CMB power spectrum_ , _JCAP_ 01 (2011) 030 [1010.3693].
  * [187] D. Babich, P. Creminelli and M. Zaldarriaga, _The Shape of non-Gaussianities_ , _JCAP_ 08 (2004) 009 [astro-ph/0405356].
  * [188] J.R. Fergusson and E.P.S. Shellard, _The shape of primordial non-Gaussianity and the CMB bispectrum_ , _Phys. Rev. D_ 80 (2009) 043510 [0812.3413].
  * [189] Planck collaboration, _Planck 2018 results. IX. Constraints on primordial non-Gaussianity_ , _Astron. Astrophys._ 641 (2020) A9 [1905.05697].
  * [190] T. Suyama and M. Yamaguchi, _Non-Gaussianity in the modulated reheating scenario_ , _Phys. Rev. D_ 77 (2008) 023505 [0709.2545].
  * [191] K. Matsumoto, _Appell and lauricella hypergeometric functions_ , in _Encyclopedia of Special Functions: The Askey-Bateman Project_ , T.H. Koornwinder and J.V. Stokman, eds., vol. 2, p. 79–100, Cambridge University Press (2020), DOI.