# Cosmological Correlators with Double Massive Exchanges: Bootstrap Equation and Phenomenology

Shuntaro Aoki, Lucas Pinol, Fumiya Sano, Masahide Yamaguchi, Yuhang Zhu

► To cite this version:

Shuntaro Aoki, Lucas Pinol, Fumiya Sano, Masahide Yamaguchi, Yuhang Zhu. Cosmological Correlators with Double Massive Exchanges: Bootstrap Equation and Phenomenology. Journal of High Energy Physics, 2024, 2024 (09), pp.176. 10.1007/JHEP09(2024)176. hal-04576598

HAL Id: hal-04576598

https://hal.science/hal-04576598v1

Submitted on 10 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Cosmological correlators with double massive exchanges: bootstrap equation and phenomenology

Shuntaro Aoki,ª Lucas Pinol,ᵇ Fumiya Sano,ᶜ,>d Masahide Yamaguchi,>d,c
and Yuhang Zhu,>d

E-mail: shuntaro1230@gmail.com, lucas.pinol@phys.ens.fr, sanof.cosmo@gmail.com, gucci@ibs.re.kr, yhzhu@ibs.re.kr

**ABSTRACT:** Using the recently developed cosmological bootstrap method, we compute the exact analytical solution for the seed integral appearing in cosmological correlators with *double massive scalar exchanges*. The result is explicit, valid in any kinematic configuration, and free from spurious divergences. It is applicable to any number of fields' species with any masses. With an appropriate choice of variables, the results contain only single-layer summations. We also propose simple approximate formulas valid in different limits, enabling direct and instantaneous evaluation. Supported by exact numerical results using CosmoFlow, we explore the phenomenology of double massive exchange diagrams. Contrary to single-exchange diagrams with ubiquitous Lorentz-covariant interactions, the size of the cubic coupling constant can be large while respecting perturbativity bounds. Because of this property, the primordial bispectrum from double-exchange diagrams can be as large as, coincidentally, current observational constraints. In addition to being sizable on equilateral configurations, we show that the primordial bispectrum exhibits a large cosmological collider signal in the squeezed limit, making the double massive exchanges interesting channels for the detection of massive primordial fields. We propose to decisively disentangle double-exchange channels from single-exchange ones with cosmological observations by exploiting the phase information of the cosmological collider signal, the inflationary flavor oscillations from multiple fields' species exchanges and the double soft limit in the primordial trispectrum.

**KEYWORDS:** de Sitter space, Space-Time Symmetries, Early Universe Particle Physics

**ARXIV EPRINT:** 2404.09547

<footnote>ªParticle Theory and Cosmology Group, Center for Theoretical Physics of the Universe, Institute for Basic Science, Daejeon, 34126, Korea</footnote>

<footnote>ᵇLaboratoire de Physique de l'École Normale Supérieure, ENS, CNRS, Université PSL, Sorbonne Université, Université Paris Cité, F-75005, Paris, France</footnote>

<footnote>ᶜDepartment of Physics, Tokyo Institute of Technology, Tokyo, 152-8551, Japan</footnote>

<footnote>ᵈCosmology, Gravity and Astroparticle Physics Group, Center for Theoretical Physics of the Universe, Institute for Basic Science, Daejeon, 34126, Korea</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Seed integral and bootstrap equations</td><td>6</td></tr><tr><td>2.1</td><td>Definition of the seed integral and typical examples</td><td>7</td></tr><tr><td>2.2</td><td>Derivation of bootstrap equations with double-exchange</td><td>9</td></tr><tr><td>3</td><td>Analytical computations for bootstrap equations</td><td>15</td></tr><tr><td>3.1</td><td>Particular solutions</td><td>16</td></tr><tr><td>3.2</td><td>Determine the coefficients</td><td>21</td></tr><tr><td>3.3</td><td>Summary of solutions</td><td>21</td></tr><tr><td>4</td><td>Continuation to all kinematic regions</td><td>23</td></tr><tr><td>4.1</td><td>Continuation and resummation</td><td>23</td></tr><tr><td>4.2</td><td>Consistency checks at different limits</td><td>27</td></tr><tr><td>5</td><td>Double massive exchange phenomenology</td><td>31</td></tr><tr><td>5.1</td><td>Equilateral value and mass dependence</td><td>36</td></tr><tr><td>5.2</td><td>Shape dependence</td><td>36</td></tr><tr><td>5.3</td><td>Cosmological collider signal</td><td>38</td></tr><tr><td>5.4</td><td>Disentangling double-exchange from single-exchange channels</td><td>43</td></tr><tr><td>6</td><td>Conclusions and outlooks</td><td>49</td></tr><tr><td>A</td><td>Mathematical formulae</td><td>51</td></tr><tr><td>B</td><td>Seed integral with single massive exchange</td><td>53</td></tr><tr><td>C</td><td>Boundary conditions of seed integrals</td><td>54</td></tr><tr><td>D</td><td>Toward triple exchange</td><td>56</td></tr></table>

## 1 Introduction

**Cosmological correlators.** Key predictions of inflationary cosmology are encoded as statistical information in the cosmological correlators, namely, the *n*-point correlation functions of ubiquitous massless fields: the curvature perturbation and the gravitons. The final values of these correlation functions at the end of inflation depend on the dynamics and interactions of quantum fluctuations in the bulk of the inflationary spacetime. Because these correlators provide the initial conditions for the subsequent evolution of the universe, observing the correlation functions in the Cosmic Microwave Background (CMB) and Large Scale Structures (LSS) offers the intriguing possibility of gaining insights into high-energy processes inaccessible through terrestrial collider experiments. The correlators involving the exchange of massive
and possibly spinning particles are of particular interest. The inherently high-energy scale of inflation creates an environment favorable for producing massive fields during this remarkable epoch. Through their interactions with the curvature perturbation, these fields generate distinct patterns in the correlation functions, encoding crucial information about their masses, spins, mixing angles, and interactions. This concept was initially introduced in pioneering works [1–3], and was subsequently elaborated and termed as the *Cosmological Collider* (CC) [4]. Numerous recent studies have utilized the CC signal to investigate models of high-energy physics [5–99] along with developing both analytical and numerical techniques aiming to deepen the understanding of underlying fundamental properties of these correlators [100–155], see also the recent first searches for CC signals in real cosmological data [156, 157].

**Analytical and numerical methods.** Conventional methods for calculating inflationary correlation functions, based on the so-called in-in (or Schwinger-Keldysh (SK)) formalism, typically involve multiple layers of nested integrals in the time domain (see [158–160] for comprehensive reviews). Remarkably, in the recent years we have seen significant developments in various analytical and numerical methods. For instance, newly developed wavefunction approaches have made the analytical properties behind correlators more transparent [110–114, 124, 133, 150, 161–169]. Indeed, evaluating directly the nested time integrals — involving massive propagators expressed as complicated special functions — can be extremely challenging. Given that the inflationary background is quasi-de Sitter (dS), by restricting the calculation of correlators to their values at the end of inflation only, we can bring a drastic simplification. This so-called “boundary perspective” was introduced in a program known as the *cosmological bootstrap* [106–108], see [128] for a review. Initially relying on the exact dS symmetries, this method has led to the full analytical form of tree-level correlators built off single-exchange channels only. Later, it has been extended to symmetry-breaking cases, such as those with non-unit sound speed [130–132], chemical potential [135], with IR effects [137], and time-dependent mass [147]. Additionally, fundamental properties including unitarity, locality, causality, and analyticity were invoked to comprehend the underlying structure of these cosmological correlators [110–116]. The Mellin transformation also proved to be powerful for understanding inflationary correlators [121–123]. The so-called partially Mellin-Barnes (MB) method was introduced over the last two years [82, 135, 138–141] and makes it easier to calculate bulk time integrals, and subsequently many exact results for inflationary correlators were obtained. In a different direction, a systematic program called the “Cosmological Flow” was proposed last year [170, 171]. It enables us to trace the time evolution of the correlators in the bulk of the inflationary evolution and is valid even in the strong quadratic mixing regime that remains so far inaccessible with traditional methods. The numerical implementation of this method is called *CosmoFlow* and a user-friendly guide with non-trivial examples is provided in [172].

**Beyond single-exchange diagrams.** After these efforts, the analytical structure of correlators associated with a single massive (spinning) field exchange is now well understood. Even in theories involving several different species of massive fluctuations with non-trivial interactions, the single-exchange channel was computed, and a new phenomenon, dubbed inflationary flavor oscillations, was discovered [77]. Nevertheless, *More is different*. Moving beyond the
single-exchange tree-level diagrams is highly challenging due to the significant increase in complexity. For example, our understanding of loop diagrams is still far from satisfactory. Recently, some progress has been made towards understanding the structure of loop diagrams involving massive fields. For example, using the techniques of spectral decomposition in dS to obtain the one-loop results with the pair of massive scalars [136], and using the MB representation to extract the non-analytical signals within loop diagrams [139, 140]. Another important but complicated case arises from including several massive propagators and more interactions in the tree-level diagrams. Those diagrams are crucial because of their rich phenomenology, whereas their analytical understanding is still limited. From the perspective of the bulk time integral, incorporating more massive fields and more interaction vertices results in exceedingly complicated integrals. The integrand comprises products of several special functions, and the number of layers of nested time integrals has increased. Thus, it seems hopeless to perform the time integration directly. To address the problem, several possible methods have been proposed, each with their own unique advantages. First, an approximate method, called the bulk-version cutting rule, has proven to be useful. This method treats one massive propagator at a time while mimicking others using some local operators in an effective way, which enables us to extract the leading-order CC signals [125]. However, sub-leading contributions, as well as the background signal that is non-negligible in kinematic configurations relevant for comparison with cosmological observations, are missed. Regarding exact calculations, the Cosmological Flow and its numerical implementation CosmoFlow are useful tools as they can give the exact result for any massive spinless exchange channel, including CC and background signals, and also including the regime of strong quadratic mixing [170–172]. As far as exact analytical methods are concerned, only one very recent work has successfully generalized the MB method to include more than just single-exchange correlators [141], although the exact results are formulated with several layers of series summation. As for the bootstrap equations, to the best of our knowledge, no result has been published so far that is associated to diagrams beyond the single exchange.

**Content of this work.** In this work, we present the first success in applying the cosmological bootstrap equations for deriving the exact solutions of correlators that involve *double* massive fields exchanges. Moreover, we perform for the first time a detailed comparison of the bootstrap result with the Cosmological Flow method — completely independent since not relying on a perturbative scheme in terms of the quadratic mixing — by utilizing CosmoFlow. To be more specific, we consider the double-exchange four- and three-point functions depicted in figure 1. In full generality, we allow those two exchanged fields to have different masses, and the obtained results are applicable to theories with various types of interactions, as long as they do not lead to secular divergences. Although the concrete calculation concerns correlators of massless scalar fluctuations, introduced as external fields, the extension to graviton correlators is straightforward.

*Bootstrap equations.* The first part of our work focuses on deriving the bootstrap equations satisfied by the *seed integral* associated with these diagrams, and obtaining analytical solutions for the resulting differential equations. More precisely, we can apply a differential operator to each massive propagator individually, resulting in two second-order differential equations
with a source term being the single-exchange channel in the absence of the corresponding propagator. This whole procedure is summarized in the schematic diagram in figure 2, and is fully explained in section 2. Here, we stand on the shoulders of giants, benefiting from the well-established result of the single-exchange diagram, particularly the exact and closed-form expression of the corresponding seed integral [138] which does not require any series summation.

**Solutions to the bootstrap equations.** We show in section 3 that we can solve these complicated bootstrap equations. The homogeneous equations can be solved using a specific type of two-variable hypergeometric function known as the Appell function $F_4$. We have observed that a particular choice of variables is quite useful, namely $r_i \equiv 2k_i/(k_1+k_2+k_3+k_4)$, where $k_i = |k_i|$. By employing these variables and choosing appropriate ansatz, we find particular solutions and also the integration constants, ultimately deriving the exact general solutions for these bootstrap equations. Notably, these solutions involve only *one* series summation which makes the analysis of their analytical structure more transparent. As is now well known, the single-exchange channel exhibits two distinct contributions characterized by different behaviors: one contribution involves imaginary powers of momentum ratios and is known as the celebrated CC signal, while the other contributions display analytic dependence on the momentum ratios and are typically regarded as the featureless background. As could be expected, the double-massive fields exchanges channels consist of three distinct terms: one arises from signals generated by both massive fields, another one from the mixing between the background and the signal originating from only one massive field and the remaining one represents purely the background. Each term has different mass dependence as will be emphasized in the main text. This iterative procedure, employing simpler diagrams as the source to bootstrap equations for seed integrals involving more massive fields exchanges, can be generalized to more complicated diagrams. In particular, we demonstrate that the results of the double-exchange diagram can be used to calculate triple-exchange correlators.

**Explicit cancellations of spurious divergences.** Unfortunately, the primary results of the four-point correlation function are not applicable to all kinematic regions, as also observed in [141]. The situation is exacerbated when considering the three-point function, typically derived by taking the soft limit of one of the external lines (e.g., $k_4 \to 0$ in figure 1). The primary result obtained from taking this limit does not converge for kinematic regions that satisfy the triangle inequality. In section 4, we will perform some transformations of the special functions in it, and carefully examine the divergent behaviour under various limits (e.g. folded limit). We observed that, by regrouping the divergences of homogeneous and particular solutions together, the whole expression achieves convergent. In brief, by using the continuation of certain special functions and appropriately organizing and treating all contributions as a whole, we are able to obtain a final result that is explicitly convergent and valid for all physical kinematic regions, rendering it perfectly usable for concrete purposes.

**Large primordial bispectrum and cosmological collider signal, inflationary flavor oscillations and the primordial trispectrum.** The second part of our work, in section 5, precisely digs into the rich phenomenology associated with double massive exchange (DE) diagrams. We
recall how to relate correlators of the massless scalar fluctuation $\phi$ to cosmological correlation functions of the primordial curvature fluctuation $\zeta$. We also propose two natural embeddings for the theory under study, with the massive scalar fluctuations $\sigma_\alpha$ identified as the massive eigenstates of the isocurvature fluctuations during inflation. We explain that, in the large class of general non-linear sigma models of inflation, one should generically expect cubic interactions involving one massless field and two massive ones, leading to DE diagrams. Moreover, from generic effective field theory perspectives, we also expect such massive fluctuations to be coupled via an interaction of this form. Importantly, the size of this cubic interaction — set by the curvature of the target space in non-linear sigma models — is not dictated by the quadratic mixing alone, contrary to the ubiquitous Lorentz-covariant one leading to single-exchange (denoted as SE in the following) diagrams. Therefore, we show that the resulting primordial bispectrum can be large, first on equilateral configurations for which we display a simple fitting formula, and also in squeezed configurations for which we derive the CC signal in a closed and explicit form for the first time. Importantly, the DE CC signal is not suppressed by additional factors of the Boltzmann factor ($e^{-\pi\mu}$ where $\mu = \sqrt{m^2/H^2 - 9/4}$ is the mass parameter of the double exchanged massive field) compared to the SE CC signal, and we show that its relative polynomial suppression ($\mu^{-2}$) can easily be compensated by a larger cubic coupling constant for masses $m \gtrsim H$ relevant for cosmological observations. We explain that the bootstrap result is exact and valid for any kinematic configuration, showcasing its usefulness to predict the full shape of the bispectrum and making contact with observations. For each of these steps, we precisely compare our analytical predictions with exact numerical calculations using CosmoFlow and not relying on a perturbative scheme for the quadratic mixing. Impressively, the two independent methods yield almost identical results for any kinematic configuration, and any values of the mass of the double exchanged field. Moreover, supported by the exact numerical result, we propose a first naive extrapolation of the bootstrap result to the strong quadratic mixing regime, highlighting the weaknesses of this procedure. We also propose ways to disentangle DE channels from SE ones from cosmological observations, beyond the simple observation that the DE signal may be observable while respecting perturbativity bounds, contrary to the SE one. First, we showcase the utility of the CC phase information to tell those channels apart, since the relation between frequency and phase are well different for masses close to the Hubble scale. Second, we explain that in the more generic situation where there exist different species of massive fields, each with their own masses, one can use the inflationary flavor oscillations to tell apart the two channels. Third, by moving to the primordial trispectrum we show a unique feature of DE diagrams that cannot be mimicked by a SE one, with CC signal oscillations transitioning from a frequency $\mu$ to $2\mu$ in the double soft limit.

**Organization of this work.** This paper is organized as follows. In section 2, we define the seed integral related to the double-exchange inflation correlators and then derive the bootstrap equations satisfied by these seed integrals. In section 3, we derive the exact solution of the bootstrap equations step by step, and the final exact expression is summarized in section 3.3. The results are extended to all physical regions in section 4, where we also discuss their behaviors under various limits. In section 5, we delve into the intriguing phenomenology underlying the double-exchange correlators. Finally, we conclude and outline
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Four-point and three-point diagrams with double massive exchange. For generality, here we allow internal massive scalar fields $\sigma_{\alpha}$ and $\sigma_{\beta}$ to have different masses.</figcaption></figure>

future directions in section 6.

**Conventions and notations.** Throughout this paper, we adopt the $(-, +, +, +)$ metric sign convention. The background spacetime is fixed as: $ds^2 = a^2(\tau)(-d\tau^2 + dx^2)$, with the scale factor $a(\tau) = -1/H\tau$ and $\tau \in (-\infty, 0)$. Here $H$ represents the Hubble parameter that is fixed as the constant. A prime on quantum average values $\langle \cdots \rangle'$ indicates the momentum $\delta$-function and the factor $(2\pi)^3$ is omitted. We frequently use the shorthand notations, for example $k_{1234} \equiv \sum_{i=1}^4 k_i$ and $p_{123} \equiv p_1 + p_2 + p_3$, with similar shorthand subscripts following the same convention. Some frequently used variables are defined as

$$u \equiv \frac{k_1}{k_{24}}, \quad v \equiv \frac{k_3}{k_{24}}, \quad r_i \equiv \frac{2k_i}{k_{1234}}, \quad \tilde{r}_i \equiv \frac{2k_i}{k_{123}}.$$

The majority of the special functions utilized in this work are referenced in the mathematical functions handbook [173]. Additionally, in appendix A, we provide the definitions and useful formulae related to these special functions. Other variables and functions will be defined in the main text as they are introduced.

## 2 Seed integral and bootstrap equations

This part serves as an introduction to the basic ingredients needed for later discussion. Experts who are familiar with these descriptions can directly skip to section 2.1.

As mentioned in the introduction, in this work, we will mainly focus on the four-point and three-point inflation correlators shown in figure 1.

The external lines here always represent the massless inflaton fluctuation $\phi$ (or, equivalently, the curvature perturbation $\zeta$), and the internal lines are associated with massive scalar fields $\sigma_{\alpha}$ ($\alpha = 1, 2, \cdots$). Distinguished by different colors for clarity, we allow the exchange of massive scalar fields with distinct masses $m_{\alpha}$, where $\alpha$ (or $\beta$) serves as the label for different species. By taking the free Hamiltonian to be quadratic and diagonal, each of the fields can be quantized separately as

$$\Phi^A(\tau, \mathbf{x}) = \int \frac{d^3\mathbf{k}}{(2\pi)^3} (\Phi_k^A(\tau) a_{\mathbf{k}}^A + \Phi_k^{A*}(\tau) a_{-\mathbf{k}}^{A\dagger}) e^{i\mathbf{k}\cdot\mathbf{x}} \quad (\text{no sum on } A), \quad (2.1)$$

where $\Phi^A$ can denote either the inflaton perturbation $\phi$ or the massive fields $\sigma_{\alpha}$ and $a_{\mathbf{k}}^{A\dagger}(a_{\mathbf{k}}^A)$ are the creation (annihilation) operators satisfying the usual commutation relations $[a_{\mathbf{k}}^A, a_{\mathbf{k}'}^{A\dagger}] =$
(2π)$^3\delta^{AB}\delta^{(3)}(\mathbf{k}-\mathbf{k}')$. Mode functions $\Phi_k^A(\tau) = \{\varphi_k, \sigma_k^\alpha\}$ are specified below,

$$ \varphi_k = \frac{H}{\sqrt{2k^3}}(1 + ik\tau)e^{-ik\tau}, \qquad (2.2) $$

$$ \sigma_k^\alpha = -i e^{-\frac{\pi}{2}\mu_\alpha + i\frac{\pi}{4}} \frac{\sqrt{\pi}}{2} H(-\tau)^{3/2} H_{i\mu\alpha}^{(1)}(-k\tau), \qquad (2.3) $$

where $\mu_\alpha = \sqrt{m_\alpha^2/H^2 - 9/4}$ and $H_\nu^{(1)}(\cdot)$ is the Hankel function of the first kind of order $\nu$. In all the following sections, we show our results for the case with massive fields in the principal series ($m_\alpha > 3H/2$) which is particularly relevant to cosmological collider physics. The results related to imaginary $\mu_\alpha$ ($m_\alpha < 3/2H$) can then be simply obtained by the substitution $\mu_\alpha \to -i\nu_\alpha$ with $\nu_\alpha = \sqrt{9/4 - m_\alpha^2/H^2}$, and we will use that for a specific example later on in section 5. To formulate the seed integral, we adopt the Schwinger-Keldysh (SK) diagrammatic conventions outlined in [160], where readers can find more details. The bulk-to-boundary propagators $K$ associated with the inflaton fluctuation $\varphi$ are

$$ K_a(k, \tau) = \frac{H^2}{2k^3} (1 - iak\tau) e^{iak\tau}, \qquad (2.4) $$

where $a = \pm 1$ is the SK vertex index, and the bulk-to-bulk propagators $D_{ab}^\alpha$ associated with massive fields $\sigma_\alpha$ are explicitly given by

$$ D_{+-}^{\alpha}(k; \tau_1, \tau_2) = \sigma_k^{\alpha}(\tau_1) \sigma_k^{\alpha*}(\tau_2) = \frac{H^2 \pi e^{-\pi \mu_{\alpha}}}{4} (\tau_1 \tau_2)^{3/2} H_{i\mu\alpha}^{(1)}(-k\tau_1) H_{-i\mu\alpha}^{(2)}(-k\tau_2), \quad (2.5) $$

$$ D_{+-}^{\alpha}(k; \tau_1, \tau_2) = (D_{+-}^{\alpha}(k; \tau_1, \tau_2))^*, \qquad (2.6) $$

$$ D_{\pm\pm}^{\alpha}(k; \tau_1, \tau_2) = D_{\mp\pm}^{\alpha}(k; \tau_1, \tau_2) \theta(\tau_1 - \tau_2) + D_{\pm\mp}^{\alpha}(k; \tau_1, \tau_2) \theta(\tau_2 - \tau_1). \qquad (2.7) $$

Here $\tau_1$ and $\tau_2$ denote the conformal time at the vertex where bulk-to-bulk propagators connect. Note that propagators $D_{\pm\pm}^\alpha$ are factorised in time and satisfy the homogeneous differential equation, whereas the time ordering in the $D_{\pm\pm}$ introduced an additional $\delta$-function source to this equation, explicitly

$$ \left[ \tau_i^2 \partial_{\tau_i}^2 - 2\tau_i \partial_{\tau_i} + k^2 \tau_i^2 + \mu_\alpha^2 + \frac{9}{4} \right] D_{\pm\mp}^\alpha (k; \tau_1, \tau_2) = 0, \qquad (2.8) $$

$$ \left[ \tau_i^2 \partial_{\tau_i}^2 - 2\tau_i \partial_{\tau_i} + k^2 \tau_i^2 + \mu_\alpha^2 + \frac{9}{4} \right] D_{\pm\pm}^\alpha (k; \tau_1, \tau_2) = \mp i H^2 \tau_1^2 \tau_2^2 \delta(\tau_1 - \tau_2), \qquad (2.9) $$

where $i = 1, 2$. These differential equations will play an important role in deriving the bootstrap equations later and one can easily check that they are indeed satisfied by (2.3).

## 2.1 Definition of the seed integral and typical examples

Based on the quantities defined in the previous sections, we now introduce a seed integral $\mathcal{I}$ for four-point inflaton correlators involving double massive field exchange:

$$ \mathcal{I}_{abc, \alpha\beta}^{p_1 p_2 p_3} = H^{-4} k_{24}^{9+p_{123}} (-iabc) \int_{-\infty}^{0} d\tau_1 d\tau_2 d\tau_3 (-\tau_1)^{p_1} (-\tau_2)^{p_2} (-\tau_3)^{p_3} \\ \times e^{iak_1\tau_1 + ibk_{24}\tau_2 + ick_3\tau_3} D_{ab}^{\alpha}(k_1; \tau_1, \tau_2) D_{bc}^{\beta}(k_3; \tau_2, \tau_3), \qquad (2.10) $$
where $p_i$ ($i = 1, 2, 3$) are constant numbers specifying types of $\varphi-\sigma$ vertices and $k_j$ ($j = 1, 2, 3, 4$) are the four external momenta of inflatons $\varphi_{k_j}$ and the prefactor $H^{-4}$ is introduced to render the seed integral dimensionless. We use the abbreviation $p_{123} \equiv p_1 + p_2 + p_3$, $k_{24} \equiv k_2 + k_4$, and assume $p_i > -5/2$ for the time-integrals to converge at late times $\tau_i \to 0$.<fnref n="1" /> We also remind $a, b, c = \pm$ correspond to the SK indices and $\alpha, \beta = 1, 2, \dots$ label the massive fields. The propagators for massive fields $D_{ab}^\alpha$ are defined in (2.5)–(2.7). The seed integral (2.10) can be diagrammatically shown on the left of figure 1, and other channels can be obtained by simply taking the permutation of the momentum variables.

To get an analytical expression of the seed integral (2.10) is one of the main purposes of this paper. We derive the differential equations (bootstrap equations) for eq. (2.10) in the next part section 2.2 and solve them in section 3. Furthermore, by taking the limit $k_4 \to 0$ in (2.10), one can relate it to the expression for the three-point function (as depicted in the schematic diagram of the right of figure 1). Taking this limit is a non-trivial task due to cancellations of apparent divergences, and therefore this will be discussed separately in section 4.

**Example.** Once the analytical expressions for the seed integral (2.10) and its various limits are obtained, we then have all the necessary information to compute two-, three-, and four-point correlators with double exchanges for various $\varphi-\sigma$ interactions. Among them, the three-point correlators with interactions corresponding to $(p_1, p_2, p_3) = (-2, -2, -2)$ are of particular phenomenological interest since they generically arise from UV motivated multi-field inflation models with curved field space metric [174, 175]. As we will show in greater details in section 5, this interaction can also lead to a bispectrum sufficiently large to be constrained by current and upcoming cosmological observations while remaining under theoretical control, thus making it an interesting channel for discovery of new physics at high energies.

Let us look more at this example (three-point correlators with $(p_1, p_2, p_3) = (-2, -2, -2)$) and how it is expressed by the seed integral (2.10). In this case, the interactions between inflaton's fluctuations $\varphi$ and massive scalars $\sigma_\alpha$ are given by [77, 175]

$$S_{\text{int},2} = \sum_{\alpha} \int d\tau d^3\mathbf{x} \ \rho_{\alpha} a^3 \sigma_{\alpha} \varphi', \qquad (2.11)$$

$$S_{\text{int},3} = \sum_{\alpha,\beta} \int d\tau d^3\mathbf{x} \ \lambda_{\alpha\beta} a^3 \sigma_{\alpha} \sigma_{\beta} \varphi', \qquad (2.12)$$

where $\rho_\alpha$ and $\lambda_{\alpha\beta}$ are coupling constants with mass dimensions 1 and 0 respectively, and the prime denotes the derivative with respect to conformal time $\tau$. Note that $a = (-H\tau)^{-1}$ in de Sitter spacetime. Then, by the SK diagrammatic rule [160], one can write down the three-point correlator of double-exchange with interactions (2.11) and (2.12) as

$$\langle \varphi_{\mathbf{k}_1} \varphi_{\mathbf{k}_2} \varphi_{\mathbf{k}_3} \rangle' = -2i \sum_{\alpha, \beta} \frac{\rho_{\alpha} \rho_{\beta} \lambda_{\alpha\beta}}{H^9} \int_{-\infty}^{0} d\tau_1 d\tau_2 d\tau_3 (-\tau_1)^{-3} (-\tau_2)^{-3} (-\tau_3)^{-3} \\ \times \sum_{\mathbf{a}, \mathbf{b}, \mathbf{c} = \pm} (\mathbf{a}\mathbf{b}\mathbf{c}) K'_{\alpha}(k_1; \tau_1) K'_{\beta}(k_2; \tau_2) K'_{\gamma}(k_3; \tau_3) D_{\mathbf{ab}}^{\alpha}(k_1; \tau_1, \tau_2) D_{\mathbf{bc}}^{\beta}(k_3; \tau_2, \tau_3) + 2 \text{ perms}, \qquad (2.13)$$

<footnote n="1">¹As long as Re[$p_i$] > -5/2, it can be analytically continued to imaginary $p_i$ in principle [138].</footnote>
where the prime on the left-hand side means that a momentum conservation factor $(2\pi)^3\delta^{(3)}(\mathbf{k}_1+\mathbf{k}_2+\mathbf{k}_3)$ is extracted, and “2 perms” represents the permutations $(k_2 \leftrightarrow k_1)$ and $(k_2 \leftrightarrow k_3)$. The propagators $K_a$ and $D_{ab}^\alpha$ are given in (2.4) and (2.5)–(2.7), respectively. With $K'_a(k, \tau) = H^2\tau/2k e^{iak\tau}$, this can be further simplified as

$$ \langle \varphi_{\mathbf{k}_1} \varphi_{\mathbf{k}_2} \varphi_{\mathbf{k}_3} \rangle' = 2i \sum_{\alpha, \beta} \frac{\rho_{\alpha} \rho_{\beta} \lambda_{\alpha\beta}}{H^3} \cdot \frac{1}{8k_1 k_2 k_3} \int_{-\infty}^{0} d\tau_1 d\tau_2 d\tau_3 (-\tau_1)^{-2} (-\tau_2)^{-2} (-\tau_3)^{-2} \\ \times \sum_{\mathbf{a}, \mathbf{b}, \mathbf{c} = \pm} (\mathbf{a}\mathbf{b}\mathbf{c}) e^{i\mathbf{k}\mathbf{a}\mathbf{1}\tau_1 + i\mathbf{b}\mathbf{k}\mathbf{2}\tau_2 + i\mathbf{c}\mathbf{k}\mathbf{3}\tau_3} D_{\mathbf{ab}}^{\alpha}(k_1; \tau_1, \tau_2) D_{\mathbf{bc}}^{\beta}(k_3; \tau_2, \tau_3) + 2 \text{ perms.} \tag{2.14} $$

Now, it is clear that this three-point correlator can be directly related to the seed integral (2.10) with

$$ \langle \varphi_{\mathbf{k}_1} \varphi_{\mathbf{k}_2} \varphi_{\mathbf{k}_3} \rangle' = - \sum_{\alpha, \beta} \rho_{\alpha} \rho_{\beta} \lambda_{\alpha\beta} \cdot \frac{H}{4k_1 k_2^4 k_3} \sum_{\mathbf{a}, \mathbf{b}, \mathbf{c} = \pm} \lim_{k_4 \to 0} \mathcal{I}_{\mathbf{abc}, \alpha\beta}^{-2-2-2} + 2 \text{ perms.} \tag{2.15} $$

Another example that will be discussed in this work is the trispectrum with four-point interaction like

$$ S_{\text{int},4} = \sum_{\alpha, \beta} \int d\tau d^3\mathbf{x} \, \tilde{\lambda}_{\alpha\beta} a^2 \sigma_{\alpha} \sigma_{\beta} (\varphi')^2, \tag{2.16} $$

where $\tilde{\lambda}_{\alpha\beta}$ is a coupling constant with mass dimension $-2$. Following the same procedure as before, the trispectrum involving two massive propagators can be expressed as

$$ \langle \varphi_{\mathbf{k}_1} \varphi_{\mathbf{k}_2} \varphi_{\mathbf{k}_3} \varphi_{\mathbf{k}_4} \rangle' = -4i \sum_{\alpha, \beta} \rho_{\alpha} \rho_{\beta} \tilde{\lambda}_{\alpha\beta} \cdot \frac{1}{16k_1 k_2 k_3 k_4} \int_{-\infty}^{0} d\tau_1 d\tau_2 d\tau_3 (-\tau_1)^{-2} (-\tau_2)^{0} (-\tau_3)^{-2} \\ \times \sum_{\mathbf{a}, \mathbf{b}, \mathbf{c} = \pm} (\mathbf{a}\mathbf{b}\mathbf{c}) e^{i\mathbf{k}\mathbf{a}\mathbf{1}\tau_1 + i\mathbf{b}\mathbf{k}\mathbf{2}\tau_2 + i\mathbf{c}\mathbf{k}\mathbf{3}\tau_3} D_{\mathbf{ab}}^{\alpha}(k_1; \tau_1, \tau_2) D_{\mathbf{bc}}^{\beta}(k_3; \tau_2, \tau_3) + 5 \text{ perms.} \tag{2.17} $$

where the total of $6=4!/4$ permutations correspond to the different ways to form unordered pairs out of the four $k_i$. Then, the trispectrum can be connected to the seed integral (2.10) through a simple relation, such as

$$ \langle \varphi_{\mathbf{k}_1} \varphi_{\mathbf{k}_2} \varphi_{\mathbf{k}_3} \varphi_{\mathbf{k}_4} \rangle' = \sum_{\alpha, \beta} \rho_{\alpha} \rho_{\beta} \tilde{\lambda}_{\alpha\beta} \cdot \frac{H^4}{4k_1 k_2 k_3 k_4 k_{24}^5} \sum_{\mathbf{a}, \mathbf{b}, \mathbf{c} = \pm} \mathcal{I}_{\mathbf{abc}, \alpha\beta}^{-2\ 0\ -2} + 5 \text{ perms.} \tag{2.18} $$

## 2.2 Derivation of bootstrap equations with double-exchange

Now, our task is to compute the seed integral (2.10). However, the nested time integral of the special function makes it difficult to perform a direct integration. Instead of doing so, here we derive differential equations that the seed integral satisfies and subsequently solve them with appropriate boundary conditions. To facilitate the later derivation, let us first change the integration variables from $\tau_i$ to $z_i$ by

$$ -k_1 \tau_1 = z_1, \quad -k_{24} \tau_2 = z_2, \quad -k_3 \tau_3 = z_3, \tag{2.19} $$
and hence $z_i \in (0, +\infty)$. Then, the propagator $D_{ab}^{\alpha}$ inside the seed integral (2.10) can be rewritten as

$$D_{ab}^{\alpha}(k_1; \tau_1, \tau_2) = D_{ab}^{\alpha}\left(k_1; -\frac{z_1}{k_1}, -\frac{z_2}{k_{24}}\right) \equiv k_1^{-3} \hat{D}_{ab}^{\alpha}(z_1, uz_2), \quad (2.20)$$

where we have defined momentum ratios as

$$u \equiv \frac{k_1}{k_{24}}, \qquad v \equiv \frac{k_3}{k_{24}}, \qquad (2.21)$$

and “hat” propagators $\hat{D}_{ab}$ by

$$\hat{D}_{-+}^{\alpha}(z_1, uz_2) = \frac{H^2 \pi e^{-\pi \mu_{\alpha}}}{4} z_1^{3/2} (uz_2)^{3/2} H_{i\mu_{\alpha}}^{(1)}(z_1) H_{-i\mu_{\alpha}}^{(2)}(uz_2), \quad (2.22)$$

$$\hat{D}_{+-}^{\alpha}(z_1, uz_2) = (\hat{D}_{-+}^{\alpha}(z_1, uz_2))^{*}, \quad (2.23)$$

$$\hat{D}_{\pm\pm}^{\alpha}(z_1, uz_2) = \theta(uz_2 - z_1) \hat{D}_{\mp\pm}^{\alpha}(z_1, uz_2) + \theta(z_1 - uz_2) \hat{D}_{\pm\mp}^{\alpha}(z_1, uz_2). \quad (2.24)$$

In the same way, for $D_{bc}^{\beta}$ (2.10), we have

$$D_{bc}^{\beta}(k_3; \tau_2, \tau_3) = D_{bc}^{\beta}\left(k_3; -\frac{z_2}{k_{24}}, -\frac{z_3}{k_3}\right) = k_3^{-3} \hat{D}_{bc}^{\beta}(vz_2, z_3). \quad (2.25)$$

With the redefinition above, the seed integral (2.10) is written as

$$\begin{aligned} \mathcal{I}_{abc, \alpha\beta}^{p_1 p_2 p_3} &= \frac{H^{-4}(-iabc)}{u^{4+p_1}v^{4+p_3}} \int_0^{\infty} dz_1 dz_2 dz_3 z_1^{p_1} z_2^{p_2} z_3^{p_3} e^{-iaz_1 - ibz_2 - icz_3} \hat{D}_{ab}^{\alpha}(z_1, uz_2) \hat{D}_{bc}^{\beta}(vz_2, z_3), \\ &\equiv \frac{1}{u^{4+p_1}v^{4+p_3}} \hat{\mathcal{I}}_{abc, \alpha\beta}^{p_1 p_2 p_3}(u, v), \end{aligned} \quad (2.26)$$

where we have introduced the hatted seed integral $\hat{\mathcal{I}}$ in the second line. Note that the seed integral depends on the two combinations of momentum, $u = k_1/k_{24}$ and $v = k_3/k_{24}$. The important observation here is that $\hat{D}_{ab}^{\alpha}(\hat{D}_{bc}^{\beta})$ depends on a specific combination $uz_2$ ($vz_2$), which allows us to derive differential equations for the hat propagators with respect to $u$ ($v$) instead of $\tau_2$ or $z_2$. For example, from eqs. (2.8) and (2.9) with change of variables (2.19), we find that $\hat{D}_{ab}^{\alpha}(z_1, uz_2)$ satisfy

$$\left[ z_2^2 \partial_{z_2}^2 - 2z_2 \partial_{z_2} + (uz_2)^2 + \mu_{\alpha}^2 + \frac{9}{4} \right] \hat{D}_{\pm\mp}^{\alpha}(z_1, uz_2) = 0, \quad (2.27)$$

$$\left[ z_2^2 \partial_{z_2}^2 - 2z_2 \partial_{z_2} + (uz_2)^2 + \mu_{\alpha}^2 + \frac{9}{4} \right] \hat{D}_{\pm\pm}^{\alpha}(z_1, uz_2) = \mp i H^2 z_1^2 (uz_2)^2 \delta(z_1 - uz_2). \quad (2.28)$$

Then, noting $z_2 \partial_{z_2} f(uz_2) = u \partial_u f(uz_2)$, one can rewrite the equations above to those with respect to $u$ that is

$$\left[ u^2 \partial_u^2 - 2u \partial_u + (uz_2)^2 + \mu_{\alpha}^2 + \frac{9}{4} \right] \hat{D}_{\pm\mp}^{\alpha}(z_1, uz_2) = 0, \quad (2.29)$$

$$\left[ u^2 \partial_u^2 - 2u \partial_u + (uz_2)^2 + \mu_{\alpha}^2 + \frac{9}{4} \right] \hat{D}_{\pm\pm}^{\alpha}(z_1, uz_2) = \mp i H^2 z_1^2 (uz_2)^2 \delta(z_1 - uz_2). \quad (2.30)$$
In the same way, for $\hat{D}_{bc}^{\beta}(vz_2, z_3)$, we obtain

$$ \left[ v^2 \partial_v^2 - 2v \partial_v + (vz_2)^2 + \mu_\beta^2 + \frac{9}{4} \right] \hat{D}_{\pm\mp}^{\beta} (vz_2, z_3) = 0, \quad (2.31) $$

$$ \left[ v^2 \partial_v^2 - 2v \partial_v + (vz_2)^2 + \mu_\beta^2 + \frac{9}{4} \right] \hat{D}_{\pm\pm}^{\beta} (vz_2, z_3) = \mp i H^2 (vz_2)^2 z_3^2 \delta(vz_2 - z_3) . \quad (2.32) $$

The seed integrals consist of massive propagators. By utilizing the differential equations mentioned above and employing integration by parts, we can derive the differential equations that the seed integral satisfies. Let us begin with the simplest case, $\hat{I}_{\pm\mp\pm}$ which involves only the non-time-ordered propagators $\hat{D}_{\pm\mp}$, where the bulk time integral at each vertex is completely factorised. By applying some derivative operators, it can be easily shown

$$ \begin{aligned} & \left[ u^2 \partial_u^2 - 2u \partial_u + \mu_\alpha^2 + \frac{9}{4} \right] \hat{I}_{\pm\mp,\alpha\beta}^{p_1 p_2 p_3} \\ &= \pm i H^{-4} \int_0^\infty dz_1 dz_2 dz_3 z_1^{p_1} z_2^{p_2} z_3^{p_3} e^{\mp i z_1 \pm i z_2 \mp i z_3} \left[ u^2 \partial_u^2 - 2u \partial_u + \mu_\alpha^2 + \frac{9}{4} \right] \hat{D}_{\pm\mp}^{\alpha} (z_1, uz_2) \hat{D}_{\mp\pm}^{\beta} (vz_2, z_3) \\ &= \mp i u^2 H^{-4} \int_0^\infty dz_1 dz_2 dz_3 z_1^{p_1} z_2^{p_2+2} z_3^{p_3} e^{\mp i z_1 \pm i z_2 \mp i z_3} \hat{D}_{\pm\mp}^{\alpha} (z_1, uz_2) \hat{D}_{\mp\pm}^{\beta} (vz_2, z_3), \end{aligned} \quad (2.33) $$

where we used homogeneous equation (2.29) from the second to the third line. The equation (2.33) is proportional to the original $\hat{I}^{p_1 p_2+2, p_3}$ where the power index $p_2$ is increased by a factor of two, and one can relate it to $\hat{I}^{p_1 p_2 p_3}$ by the following formula

$$ \int_0^\infty dz \, z^{p+2} e^{-i\alpha z} f(uz)g(vz) = -(u\partial_u + v\partial_v + p+2)(u\partial_u + v\partial_v + p+1) \int_0^\infty dz \, z^p e^{-i\alpha z} f(uz)g(vz), \quad (2.34) $$

for well-behaved functions $f(uz)$ and $g(vz)$, and with $a = \pm 1$.<fnref n="2" /> It turns out that (2.33) can be written as

$$ \left[ u^2 \partial_u^2 - 2u \partial_u + \mu_\alpha^2 + \frac{9}{4} \right] \hat{I}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3} = u^2 (u \partial_u + v \partial_v + p_2 + 2) (u \partial_u + v \partial_v + p_2 + 1) \hat{I}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3} . \quad (2.35) $$

$$ \begin{aligned} 0 &= \int_0^\infty dz \, \partial_z [z^{p+1} e^{-i\alpha z} f(uz)g(vz)] \\ &= \int_0^\infty dz \, z^p e^{-i\alpha z} [p + 1 - i\alpha z + u\partial_u + v\partial_v] f(uz)g(vz), \end{aligned} $$

by $z\partial_z f(uz) = u\partial_u f(uz)$ and so on. Therefore, we have

$$ \int_0^\infty dz \, z^{p+1} e^{-i\alpha z} f(uz)g(vz) = -i\alpha (u\partial_u + v\partial_v + p+1) \int_0^\infty dz \, z^p e^{-i\alpha z} f(uz)g(vz). $$

Placing $p \to p+1$ in the above equation, we arrive at

$$ \begin{aligned} & \int_0^\infty dz \, z^{p+2} e^{-i\alpha z} f(uz)g(vz) \\ &= -i\alpha (u\partial_u + v\partial_v + p+2) \int_0^\infty dz \, z^{p+1} e^{-i\alpha z} f(uz)g(vz) \\ &= -(u\partial_u + v\partial_v + p+2) (u\partial_u + v\partial_v + p+1) \int_0^\infty dz \, z^p e^{-i\alpha z} f(uz)g(vz). \end{aligned} $$

<footnote n="2">²The proof is as follows. First, we note</footnote>
By moving the right-hand side to the left, we arrive at the homogeneous differential equations satisfied by $\mathcal{I}_{\pm\mp\pm}$, which are free from source terms,

$$ \mathcal{D}_u^\alpha \hat{\mathcal{I}}_{\pm\mp\pm,\alpha\beta}^{p_1 p_2 p_3}(u,v) = 0, \quad (2.36) $$

where the differential operator $\mathcal{D}$ is defined as

$$ \begin{aligned} \mathcal{D}_u^\alpha &\equiv u^2 \partial_u^2 - 2u \partial_u + \mu_\alpha^2 + \frac{9}{4} - u^2 (u \partial_u + v \partial_v + p_2 + 2) (u \partial_u + v \partial_v + p_2 + 1) \\ &= u^2 (1 - u^2) \partial_u^2 - u^2 v^2 \partial_v^2 - 2u [1 + u^2 (p_2 + 2)] \partial_u - 2u^2 v (p_2 + 2) \partial_v \\ &\quad - 2u^3 v \partial_u \partial_v + \left(\mu_\alpha^2 + \frac{9}{4}\right) - u^2 (p_2 + 2) (p_2 + 1), \end{aligned} \quad (2.37) $$

giving the bootstrap equation for $\hat{\mathcal{I}}_{\pm\mp\pm}$ with respect to $u$. In the same way, one can derive the one with respect to $v$,

$$ \mathcal{D}_v^\beta \hat{\mathcal{I}}_{\pm\mp\pm,\alpha\beta}^{p_1 p_2 p_3}(u,v) = 0, \quad (2.38) $$

where the operator $\mathcal{D}_v^\beta$ is obtained by replacing $u \leftrightarrow v$ and $\alpha \leftrightarrow \beta$ in (2.37). In summary, we find that the seed integral $\hat{\mathcal{I}}_{\pm\mp\pm}$ simultaneously satisfies the homogeneous partial differential equations (2.36) and (2.38).

The similar procedure is applicable to the other seed integrals. However, a distinction arises when the derivative operators are applied to the time-ordered propagators $\hat{D}_{\pm\pm}$, where the $\delta$-function on the right-hand side of eq. (2.9) introduces additional source terms. For instance, in the case of $\hat{\mathcal{I}}_{\pm\pm\mp}$, we obtain:

$$ \begin{aligned} & \left[ u^2 \partial_u^2 - 2u \partial_u + \mu_\alpha^2 + \frac{9}{4} \right] \hat{\mathcal{I}}_{\pm\mp\pm,\alpha\beta}^{p_1 p_2 p_3} \\ &= u^2 (u \partial_u + v \partial_v + p_2 + 2) (u \partial_u + v \partial_v + p_2 + 1) \hat{\mathcal{I}}_{\pm\mp\pm,\alpha\beta}^{p_1 p_2 p_3} \\ & \quad + H^{-2} u^2 \int_0^\infty dz_1 dz_2 dz_3 z_1^{p_1+2} z_2^{p_2+2} z_3^{p_3} e^{\mp i z_1 \mp i z_2 \pm i z_3} \delta(z_1 - uz_2) \hat{D}_{\pm\mp}^\beta(vz_2, z_3), \end{aligned} \quad (2.39) $$

where the second line appears due to the $\delta$-function term in (2.30). Performing the $z_1$-integral and changing a variable by $\tilde{z}_2 = (u+1)z_2$, this source can be simplified as

$$ \frac{H^{-2} u^{p_1+4}}{(u+1)^{p_1+5}} \int_0^\infty d\tilde{z}_2 dz_3 \, \tilde{z}_2^{p_1+4} z_3^{p_3} e^{\mp i \tilde{z}_2 \pm i z_3} \hat{D}_{\pm\mp}^\beta \left( \frac{v}{u+1} \tilde{z}_2, z_3 \right). \quad (2.40) $$

Remarkably, this source term now contains only one massive propagator and two time integrals, which is actually associated to the seed integral for single massive exchange with one linear mixing vertex $\mathcal{I}_{ab,\alpha}^{p_1 p_2}$ and its exact closed analytical expression has already been investigated in details in the previous work [138] (also see appendix B for a summary of the results). To be more explicit, the seed integral with single-exchange that we will frequently utilize in this work is

$$ \mathcal{I}_{ab,\alpha}^{p_1 p_2} = H^{-2} k_3^{5+p_{12}} (-ab) \int_{-\infty}^{0} d\tau_1 d\tau_2 (-\tau_1)^{p_1} (-\tau_2)^{p_2} e^{iak_{124}\tau_1 + ibk_3\tau_2} D_{ab}^{\alpha}(k_3; \tau_1, \tau_2), \quad (2.41) $$
where the notation used here is the same as that we defined in the double-exchange seed integral.³

It turns out that the seed integral $\mathcal{I}_{ab,\alpha}^{p_1 p_2}(R)$ depends on a specific momentum ratio $R = 2k_3/k_{1234}$. By using the single-exchange seed integral $\mathcal{I}_{ab,\alpha}^{p_1 p_2}(R)$, eq. (2.40) can be expressed as (see eq. (B.4) for detailed derivation)

$$ \frac{u^{p_1+4}}{v^{p_1+5}} \mathcal{I}_{\pm\mp,\beta}^{p_1+4,p_3} \left( R = \frac{2v}{1+u+v} \right). \qquad (2.43) $$

Therefore, we obtain the bootstrap equation for $\hat{\mathcal{I}}_{\pm\pm\mp}$ with respect to $u$,

$$ \mathcal{D}_u^\alpha \hat{\mathcal{I}}_{\pm\pm\mp,\alpha\beta}^{p_1 p_2 p_3}(u,v) = \frac{u^{p_1+4}}{v^{p_1+5}} \mathcal{I}_{\pm\mp,\beta}^{p_1+4,p_3} \left( R = \frac{2v}{1+u+v} \right), \qquad (2.44) $$

which is an inhomogeneous differential equation. In addition, $\hat{\mathcal{I}}_{\pm\pm\mp}$ should also satisfy another differential equation when the operator $\mathcal{D}_v^\beta$ is applied. One difference is that $\mathcal{D}_v^\beta$ will act on the non-time ordering propagator $D_{\pm\mp}^\beta$ which has no source term contribution. Consequently, the differential equation related to $\mathcal{D}_v^\beta$ should be homogeneous. Following the same procedure, we can derive all bootstrap equations, which are summarized below:

The summary of bootstrap equations of different seed integrals.

* The fully factorised seed integral $\hat{\mathcal{I}}_{\pm\pm\mp}$:

$$ \mathcal{D}_u^\alpha \hat{\mathcal{I}}_{\pm\pm\pm,\alpha\beta}^{p_1 p_2 p_3}(u,v) = 0, \qquad (2.45) $$

$$ \mathcal{D}_v^\beta \hat{\mathcal{I}}_{\pm\pm\pm,\alpha\beta}^{p_1 p_2 p_3}(u,v) = 0. \qquad (2.46) $$

* The partially factorised partially nested seed integral $\hat{\mathcal{I}}_{\pm\pm\mp}$:

$$ \mathcal{D}_u^\alpha \hat{\mathcal{I}}_{\pm\pm\mp,\alpha\beta}^{p_1 p_2 p_3}(u,v) = \frac{u^{p_1+4}}{v^{p_1+5}} \mathcal{I}_{\pm\mp,\beta}^{p_1+4,p_3} \left( R = \frac{2v}{1+u+v} \right), \qquad (2.47) $$

$$ \mathcal{D}_v^\beta \hat{\mathcal{I}}_{\pm\pm\mp,\alpha\beta}^{p_1 p_2 p_3}(u,v) = 0. \qquad (2.48) $$

³Due to a non-linearly realized symmetry, whenever there exists a quadratic interaction of the form (2.11), it must come also with the following cubic interaction, denoted as the “Lorentz-covariant” combination

$$ S_{\text{int}} = - \sum_{\alpha} \int d\tau d^3\mathbf{x} \sqrt{-g} \frac{\rho_{\alpha}}{2f_{\pi}^2} g^{\mu\nu} \partial_{\mu} \varphi \partial_{\nu} \varphi \sigma_{\alpha}, \qquad (2.42) $$

where $f_\pi$ is an energy scale denoting the normalisation of the primordial curvature fluctuation, see section 5 for more details. With the insertion of one massive propagator and a quadratic vertex, this cubic interaction gives rise to a single-exchange bispectrum diagram related to the seed integral (2.41) via

$$ \langle \varphi_{\mathbf{k}_1} \varphi_{\mathbf{k}_2} \varphi_{\mathbf{k}_3} \rangle' = - \sum_{\alpha} \frac{\rho_{\alpha}^2 H^3}{f_{\pi}^2} \cdot \frac{1}{16 (k_1 k_2 k_3)^2} \cdot \sum_{\mathbf{a}, \mathbf{b}=\pm} \lim_{k_4 \to 0} \left[ \frac{k_{12}^2 - k_3^2}{k_3^2} \mathcal{I}_{ab}^{0-2} + i\mathbf{a} \frac{k_{12} (k_3^2 - k_1^2 - k_2^2)}{k_1 k_2 k_3} \mathcal{I}_{ab}^{-1-2} + \frac{k_3^2 - k_1^2 - k_2^2}{k_1 k_2} \mathcal{I}_{ab}^{-2-2} \right] + 2 \text{ perms.} $$
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The schematic diagram illustrates the double-exchange bootstrap equations. When acting on it with the differential operator $\mathcal{D}_u^\alpha (\mathcal{D}_v^\beta)$ associated with the massive fields $\sigma_\alpha (\sigma_\beta)$, the double-exchange diagram reduces to a single-exchange diagram without the corresponding massive propagator.</figcaption></figure>

• The fully nested seed integral $\hat{\mathcal{I}}_{\pm\pm}^A$:

$$ \mathcal{D}_u^\alpha \hat{\mathcal{I}}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3}(u,v) = \frac{u^{p_1+4}}{v^{p_1+5}} \mathcal{I}_{\pm\pm,\beta}^{p_12+4,p_3} \left( R = \frac{2v}{1+u+v} \right), \quad (2.49) $$

$$ \mathcal{D}_v^\beta \hat{\mathcal{I}}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3}(u,v) = \frac{v^{p_3+4}}{u^{p_2+5}} \mathcal{I}_{\pm\pm,\alpha}^{p_23+4,p_1} \left( R = \frac{2u}{1+u+v} \right). \quad (2.50) $$

The differential operator $\mathcal{D}_u^\alpha$ is defined in (2.37). The explicit expressions of $\mathcal{I}_{ab,\alpha}^{p_1 p_2}(R)$ can be found in (B.5) and (B.6).

Here “factorised” or “nested” is associated with the time integral. For the other partially factorised partially nested seed integral $\hat{\mathcal{I}}_{\mp\pm,\alpha\beta}(u, v)$, the expression can be easily obtained through replacing $u \leftrightarrow v$, $\alpha \leftrightarrow \beta$ and $p_1 \leftrightarrow p_3$ in the result of $\hat{\mathcal{I}}_{\pm\pm\mp,\alpha\beta}$, (2.47) and (2.48).

In figure 2, we show the schematic diagram of these differential equations. Once applying the differential operators $\mathcal{D}_u (\mathcal{D}_v)$, the double-exchange diagram is reduced to the single-exchange one for which exact closed results are already known from the previous work [138]. To be more specific, when the two SK indices of one massive propagator are opposite, applying the corresponding differential operator causes the right-hand side of figure 2 to vanish, resulting 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 single-exchange diagram.

Alternatively, as mentioned in the appendix of [131], one can apply another differential operator to both sides of figure 2 to further reduce another massive propagator. In doing so, the source becomes the simple contact one, but the trade-off is that the differential equations become of order four, and we noticed that it is then difficult to solve them. Instead, in this work, we stand on the shoulders of giants, leveraging the known results of the single-exchange diagram to derive those of the double-exchange one. By following a similar procedure, we can

<footnote n="4">⁴Similar differential equations for the fully nested integral were also derived in the appendix of [131], for a double-exchange bispectrum with a specific coupling.</footnote>
apply analogous differential operators for more complicated diagrams. In the double-exchange trispectrum depicted in figure 1 which we calculate in this paper, the interactions have two quadratic mixings along with one quartic coupling of the schematic form $\varphi^2\sigma^2$. Additionally, our method could be used to calculate another type of double-exchange four-point correlator, that generated by two cubic interactions, one $\varphi^2\sigma$ and one $\varphi\sigma^2$, and one quadratic mixing. Another important example is the seed integral for triple-exchange diagrams, with a four-point correlator generated by a quartic interaction $\varphi\sigma^3$ and three quadratic mixings. In appendix D, we present the first steps of this calculation, leaving the detailed work for a future publication.

Before moving to the solutions of those bootstrap equations, let us first address the limitations related to the procedure employed here. Indeed, although our analytical results will represent the exact solutions to the bootstrap equations we just derived, those equations were found using working hypotheses. The most crucial assumption is that the free Hamiltonian is diagonal, i.e. there is no linear mode mixing and, rather, quadratic mixings are treated as vertices of the perturbation theory. If we relax this assumption, we would need to consider a coupled system of quantum operators $\hat{\Phi}_k^A(\tau) = \sum_a \Phi_{a,k}^A(\tau) \hat{a}_k^a + \text{h.c.}$ with independent (creation and) annihilation ($\hat{a}_k^{\dagger a}$ and) $\hat{a}_k^a$ operators and defining a matrix of mode functions $\Phi_{a,k}^A(\tau)$ and related propagators. Our assumption forces us to consider the quadratic mixings of the form $\Omega_{AB}\Phi^{A\dagger}\Phi^B$, with $\Omega_{AB}$ an anti-symmetrix matrix, as quadratic vertices, therefore assuming their effect is negligible compared to the kinetic and mass terms (we can without loss of generality assume that the mass matrix is diagonal in the basis of these $\Phi^A$ [77]). This assumption imposes upper bounds on the quadratic mixings that we consider explicitly in this work, $\rho_\alpha \lesssim H$. By consistency, it also justifies our neglect of the other quadratic mixings between the massive fields $\sigma_\alpha$ identified in [175], as those would result in higher-order corrections suppressed by additional powers of small coupling constants. Having delineated the regime of validity of the analytical results will be important when we will dig into the phenomenology of double-exchange diagrams and compare predictions of the theory to observational constrains. It will also be important to compare our analytical results to an independent method not relying on this perturbative scheme for the quadratic mixings, as provided by the Cosmological Flow framework [170, 171] and its numerical implementation CosmoFlow [172]. These endeavours are undertaken in section 5 in great detail. We finish this section by mentioning another hypothesis leaving room for future works: we assumed that the inflaton fluctuations $\varphi$ have unit speed of sound, while be it in general single-field inflation or in the framework of the effective field theory of inflationary fluctuations, $c_s^2 \neq 1$ is ubiquitous. In particular, it was recently understood that a small speed of sound could result in interesting phenomenology in the soft limits of single-exchange diagrams, denoted as “low-speed collider” [90, 132], and it would be interesting to extend our results to this intriguing scenario.

### 3 Analytical computations for bootstrap equations

In this section, we will solve all the bootstrap equations derived in section 2. Readers who are not interested in technical details can directly skip to section 3.3 for the summary of final results.
Let us start from $\hat{\mathcal{I}}_{\pm\mp\pm}$, which satisfy the homogeneous partial differential equations (2.45) and (2.46), or explicitly,

$$ \left\{ \begin{aligned} & u^2 (1 - u^2) \partial_u^2 - u^2 v^2 \partial_v^2 - 2u [1 + u^2 (p_2 + 2)] \partial_u - 2u^2 v (p_2 + 2) \partial_v \\ & \quad - 2u^3 v \partial_u \partial_v + \left( \mu_\alpha^2 + \frac{9}{4} \right) - u^2 (p_2 + 2) (p_2 + 1) \end{aligned} \right\} \hat{\mathcal{I}}_{\pm\mp\pm,\alpha\beta}^{p_1 p_2 p_3} = 0, \quad (3.1) $$

$$ \left\{ \begin{aligned} & v^2 (1 - v^2) \partial_v^2 - v^2 u^2 \partial_u^2 - 2v [1 + v^2 (p_2 + 2)] \partial_v - 2v^2 u (p_2 + 2) \partial_u \\ & \quad - 2v^3 u \partial_v \partial_u + \left( \mu_\beta^2 + \frac{9}{4} \right) - v^2 (p_2 + 2) (p_2 + 1) \end{aligned} \right\} \hat{\mathcal{I}}_{\pm\mp\pm,\alpha\beta}^{p_1 p_2 p_3} = 0. \quad (3.2) $$

We find that this system of partial differential equations can be solved analytically by

$$ \hat{\mathcal{I}}_{\pm\mp\pm,\alpha\beta}^{p_1 p_2 p_3} = \sum_{a,b=\pm} c_{\pm\mp,ab} u^{\frac{3}{2}-i\alpha\mu_\alpha} v^{\frac{3}{2}-ib\mu_\beta} \mathcal{F}_4 \left[ \frac{4+p_2-i(a\mu_\alpha+b\mu_\beta)}{2} \frac{5+p_2-i(a\mu_\alpha+b\mu_\beta)}{2} \middle| u^2, v^2 \right], \quad (3.3) $$

where $c_{\pm\mp,ab}$ with $a, b = \pm$ are four undetermined coefficients, which will be fixed later by imposing proper boundary conditions. $\mathcal{F}_4[\cdots]$ is the dressed Appell series defined in (A.9). Note that $\mathcal{F}_4[\cdots]$ with this definition is only convergent for $u + v < 1$, and we will discuss in detail the continuation to all physical kinematic regions in the next section.

The other seed integrals satisfy inhomogeneous equations, and therefore, the general solutions can be expressed as the sum of homogeneous solutions, as shown in (3.3), and particular solutions, which will be discussed in the following subsections.

## 3.1 Particular solutions

Here we aim to find particular solutions for the inhomogeneous bootstrap equations represented by (2.47)–(2.50). These particular solutions are denoted as $\hat{\mathcal{P}}_{\pm\mp,\alpha\beta}^{p_1 p_2 p_3}$ and $\hat{\mathcal{P}}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3}$, respectively.

**The particular solution** $\hat{\mathcal{P}}_{\pm\mp\pm}$. This solution contributes to the partially factorised partially nested seed integral and must satisfy a system of specific equations

$$ \mathcal{D}_u^\alpha \hat{\mathcal{P}}_{\pm\mp,\alpha\beta}^{p_1 p_2 p_3} (u, v) = \frac{u^{p_1+4}}{v^{p_1+5}} \mathcal{I}_{\pm\pm,\beta}^{p_1+4, p_3} \left( R = \frac{2v}{1+u+v} \right), \quad (3.4) $$

$$ \mathcal{D}_v^\beta \hat{\mathcal{P}}_{\pm\mp,\alpha\beta}^{p_1 p_2 p_3} (u, v) = 0. \quad (3.5) $$

To find solutions, we begin by rescaling

$$ \hat{\mathcal{P}}_{\pm\pm\mp} = \frac{u^{p_1+4}}{v^{p_1+5}} \tilde{\mathcal{P}}_{\pm\pm\mp}, \quad (3.6) $$

and change variables by $r_i = 2k_i/k_{1234}$, then

$$ \frac{2u}{1+u+v} = \frac{2k_1}{k_{1234}} \equiv r_1, \quad \frac{2v}{1+u+v} = \frac{2k_3}{k_{1234}} \equiv r_3, \quad (3.7) $$
such that equations (3.4) and (3.5) become

$$ \left[ r_1^2(1-r_1)\partial_{r_1}^2 - r_1^2 r_3 \partial_{r_1} \partial_{r_3} - ((4+p_1)r_1 - 2(3+p_1))r_1 \partial_{r_1} \right. \\ \left. - (3+p_1)r_1 r_3 \partial_{r_3} + \mu_\alpha^2 + \left(p_1 + \frac{5}{2}\right)^2 \right] \tilde{\mathcal{P}}_{\pm\mp,\alpha\beta}^{p_1 p_2 p_3}(r_1, r_3) = \mathcal{I}_{\pm\mp,\beta}^{p_1 2+4, p_3}(r_3), \quad (3.8) $$

$$ \left[ r_3^2(1-r_3)\partial_{r_3}^2 - r_1 r_3^2 \partial_{r_1} \partial_{r_3} + ((5+p_{12})r_3 - 2(6+p_{12}))r_3 \partial_{r_3} \right. \\ \left. + (6+p_{12})r_1 r_3 \partial_{r_1} + \mu_\beta^2 + \left(p_{12} + \frac{13}{2}\right)^2 \right] \tilde{\mathcal{P}}_{\pm\mp,\alpha\beta}^{p_1 p_2 p_3}(r_1, r_3) = 0. \quad (3.9) $$

The differential operator appears asymmetric at the moment due to the unbalanced rescaling of $\tilde{\mathcal{P}}$ as defined in (3.6). Referring to (B.5), the right-hand side of (3.8) originates from the factorised component of the single-exchange seed integral and is explicitly provided by

$$ \mathcal{I}_{\pm\mp,\beta}^{p_1 2+4, p_3}(r_3) = \sum_{\mathbf{a}=\pm} \frac{e^{\mp i \frac{\pi}{2} (p_{12}-p_3)}}{2^{10+p_{123}-2i\mathbf{a}\mu_{\beta}} \pi^{\frac{1}{2}}} \Gamma\left[ \frac{5}{2} + p_3 + i\mu_{\beta}, \frac{5}{2} + p_3 - i\mu_{\beta}, \frac{13}{2} + p_{12} - i\mathbf{a}\mu_{\beta}, i\mathbf{a}\mu_{\beta} \middle| 3 + p_3 \right] \\ \times r_3^{\frac{13}{2}+p_{12}-i\mathbf{a}\mu_{\beta}} {}_2F_1\left[ \frac{1}{2} - i\mathbf{a}\mu_{\beta}, \frac{13}{2} + p_{12} - i\mathbf{a}\mu_{\beta} \middle| r_3 \right]. \quad (3.10) $$

The source term exhibits non-analytic behavior in $r_3$ (like a $r_3^{i\mathbf{a}\mu_{\beta}}$ term). Motivated by this observation and the series expansion of the hypergeometric function (A.3), we propose the following ansatz for $\tilde{\mathcal{P}}$,

$$ \tilde{\mathcal{P}}_{\pm\mp,\alpha\beta}^{p_1 p_2 p_3}(r_1, r_3) = \frac{e^{\mp i \frac{\pi}{2} (p_{12}-p_3)}}{2^{10+p_{123}} \pi^{\frac{1}{2}}} \Gamma\left[\frac{5}{2} + p_3 + i\mu_{\beta}, \frac{5}{2} + p_3 - i\mu_{\beta}\right] \sum_{\mathbf{a}=\pm} \sum_{m,n=0}^{\infty} \mathcal{A}_{m,n}^{\pm\pm\mp(\mathbf{a})} r_1^m r_3^{n+p_{12}+\frac{13}{2}+i\mathbf{a}\mu_{\beta}}. \quad (3.11) $$

We usually denote the oscillatory terms as CC signals, which can help identify the mass of exchanging particles based on the frequency of oscillation patterns. On the other hand, the featureless analytic terms are commonly termed as the background, upon which the signals are situated, and its effect can often be mimicked by inflaton self-interactions. Here we note that this ansatz is designed to be analytical along the $r_1$ direction but exhibit non-analytical oscillatory behavior along the $r_3$ direction as stated above. This characteristic suggests that it represents the contribution that arises from the mixing between the background and the signal. The coefficients $\mathcal{A}_{m,n}^{\pm\pm\mp(\mathbf{a})}$ are determined by the differential equations (3.8) and (3.9). Upon substituting the ansatz, we obtain the following recurrence relations,

$$ \mathcal{A}_{m+1,n}^{\pm\pm\mp(\mathbf{a})} = \frac{(3+p_1+m)(m+n+p_{12}+13/2+i\mathbf{a}\mu_{\beta})}{(m+p_1+7/2)^2 + \mu_{\alpha}^2} \mathcal{A}_{m,n}^{\pm\pm\mp(\mathbf{a})}, \quad (3.12) $$

$$ \mathcal{A}_{0,n}^{\pm\pm\mp(\mathbf{a})} = \frac{i\pi^{1/2} \operatorname{csch}(\mathbf{a}\pi\mu_{\beta})}{\mu_{\alpha}^2 + (p_1 + 5/2)^2} \Gamma\left[\frac{1}{2} + n + i\mathbf{a}\mu_{\beta}, \frac{13}{2} + p_{12} + n + i\mathbf{a}\mu_{\beta} \middle| 1 + n, \quad 1 + n + 2i\mathbf{a}\mu_{\beta}\right], \quad (3.13) $$

$$ \mathcal{A}_{m,n+1}^{\pm\pm\mp(\mathbf{a})} = \frac{(1/2 + n + i\mathbf{a}\mu_{\beta})(m + n + p_{12} + 13/2 + i\mathbf{a}\mu_{\beta})}{(n+1)(n+1 + 2i\mathbf{a}\mu_{\beta})} \mathcal{A}_{m,n}^{\pm\pm\mp(\mathbf{a})}, \quad (3.14) $$
where the first two are obtained from (3.8) while the last one is from (3.9). These recurrence relations can be solved<fnref n="5" />

$$A_{m,n}^{\pm\pm\mp(a)} = \frac{i\pi^{\frac{1}{2}} \operatorname{csch}(a\pi\mu_{\beta})}{\mu_{\alpha}^{2} + (m + p_{1} + 5/2)^{2}} \frac{(3 + p_{1})_{m}}{\left(\frac{5}{2} + p_{1} - i\mu_{\alpha}\right)_{m} \left(\frac{5}{2} + p_{1} + i\mu_{\alpha}\right)_{m}} \times \Gamma\left[\begin{array}{c} \frac{1}{2} + n + i a \mu_{\beta}, \frac{13}{2} + m + n + p_{12} + i a \mu_{\beta} \\ 1 + n, \quad 1 + n + 2 i a \mu_{\beta} \end{array}\right], \quad (3.15)$$

where $(a)_m \equiv \Gamma(a+m)/\Gamma(a)$ is the Pochhammer symbol. Upon inserting it back into (3.11), we derive the particular solution. It is noteworthy that one of the infinite summations, for instance, the n-summation, can be explicitly performed. After performing the summation over n, we obtain

$$\begin{aligned} \hat{\mathcal{P}}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3} = & \sum_{a=\pm} \frac{u^{p_1+4}}{v^{p_1+5}} \cdot \frac{i e^{\mp i \frac{\pi}{2} (p_1-p_3)}}{2^{10+p_1+5}} r_3^{p_1+\frac{13}{2}} \times \Gamma\left[\begin{array}{l} \frac{5}{2}+p_1+i\mu_\alpha, \frac{5}{2}+p_1-i\mu_\alpha, \frac{5}{2}+p_3+i\mu_\beta, \frac{5}{2}+p_3-i\mu_\beta \\ 3+p_1, 3+p_3 \end{array}\right] \\ & \times \sum_{m=0}^{\infty} \operatorname{csch}(\pi a \mu_\beta) r_1^m r_3^{i a \mu_\beta} \Gamma\left[\begin{array}{l} 3+m+p_1 \\ \frac{7}{2}+m+p_1-i\mu_\alpha, \frac{7}{2}+m+p_1+i\mu_\alpha \end{array}\right] {}_2F_1\left[\begin{array}{l} \frac{1}{2}+i a \mu_\beta, \frac{13}{2}+m+p_1+2 i a \mu_\beta \\ 1+2 i a \mu_\beta \end{array}\middle| r_3\right]. \end{aligned} \quad (3.16)$$

The particular solution $\hat{\mathcal{P}}_{\pm\pm,\alpha\beta}$. Lastly, we examine the particular solution for the fully nested seed integral $\hat{\mathcal{P}}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3}$, which satisfies the system of differential equations,

$$D_u^\alpha \hat{\mathcal{P}}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3} (u, v) = \frac{u^{p_1+4}}{v^{p_1+5}} \mathcal{I}_{\pm\pm,\beta}^{p_1+4, p_3} \left( R = \frac{2v}{1+u+v} \right), \quad (3.17)$$

$$D_v^\beta \hat{\mathcal{P}}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3} (u, v) = \frac{v^{p_3+4}}{u^{p_3+5}} \mathcal{I}_{\pm\pm,\alpha}^{p_3+4, p_1} \left( R = \frac{2u}{1+u+v} \right). \quad (3.18)$$

Similarly to the previous discussion, it is advantageous to perform rescaling in the following way

$$\hat{\mathcal{P}}_{\pm\pm,\pm} = \frac{u^{p_1+4}}{v^{p_1+5}} \tilde{\mathcal{P}}_{\pm\pm,\pm}, \quad (3.19)$$

and changing variables by (3.7). Then, the bootstrap equations are transformed into

$$\left[ r_1^2(1-r_1)\partial_{r_1}^2 - r_1^2 r_3 \partial_{r_1} \partial_{r_3} - ((4+p_1)r_1 - 2(3+p_1))r_1 \partial_{r_1} - (3+p_1)r_1 r_3 \partial_{r_3} + \mu_\alpha^2 + \left(p_1 + \frac{5}{2}\right)^2 \right] \tilde{\mathcal{P}}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3}(r_1, r_3) = \mathcal{I}_{\pm\pm,\beta}^{p_1+4, p_3}(r_3), \quad (3.20)$$

$$\left[ r_3^2(1-r_3)\partial_{r_3}^2 - r_1 r_3^2 \partial_{r_1} \partial_{r_3} + ((5+p_{12})r_3 - 2(6+p_{12}))r_3 \partial_{r_3} + (6+p_{12})r_1 r_3 \partial_{r_1} + \mu_\beta^2 + \left(p_{12} + \frac{13}{2}\right)^2 \right] \tilde{\mathcal{P}}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3}(r_1, r_3) = \left(\frac{r_3}{r_1}\right)^{p_{123}+9} \mathcal{I}_{\pm\pm,\alpha}^{p_3+4, p_1}(r_1), \quad (3.21)$$

<footnote n="5">⁵The solution is completely fixed by (3.12) and (3.13). One can explicitly check that the remaining relation (3.14) is automatically satisfied by the solution (3.15).</footnote>
where the source term originates from the nested part of the single-exchange diagram, and its exact expression can be found in (B.6),

$$ \mathcal{I}_{\pm\pm,\beta}^{p_{12}+4,p_3}(r_3) = \frac{\pm i e^{\mp i \frac{\pi}{2} p_{123}}}{2^{10+p_{123}} \pi^{\frac{1}{2}}} \Gamma\left[\frac{5}{2}+p_3+i\mu_\beta, \frac{5}{2}+p_3-i\mu_\beta\right] \times \left\{ \mp 2i \pi^{\frac{1}{2}} r_3^{p_{123}+9} {}_3F_2\left[\frac{1}{2}+p_3-i\mu_\beta, \frac{7}{2}+p_3+i\mu_\beta \middle| r_3\right] \right. \\ \left. + \sum_{\mathbf{a}=\pm} e^{\mp \pi a \mu_\beta} 2^{2ia\mu_\beta} \Gamma\left[\frac{13}{2}+p_{12}-i a \mu_\beta, i a \mu_\beta\right] r_3^{\frac{13}{2}+p_{12}-i a \mu_\beta} {}_2F_1\left[\frac{1}{2}-i a \mu_\beta, \frac{13}{2}+p_{12}-i a \mu_\beta \middle| r_3\right] \right\}, \quad (3.22) $$

where $p_\beta$ is the dressed hypergeometric function defined in (A.4), and $\mathcal{I}_{\pm\pm,\alpha}^{p_{23}+4,p_1}(r_1)$ is given by the same form as (3.22) with replacement $r_1 \leftrightarrow r_3$, $p_1 \leftrightarrow p_3$, and $\alpha \leftrightarrow \beta$. Inspired by the structure of the source term and the series expansions (A.3)-(A.4), we choose the following ansatz for the particular solutions:

$$ \tilde{\mathcal{P}}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3}(r_1, r_3) = \frac{\pm i e^{\mp i \frac{\pi}{2} p_{123}} r_3^{p_{123}}}{2^{10+p_{123}} \pi^{\frac{1}{2}}} \sum_{m,n=0}^{\infty} \left[ C_{m,n}^{\pm\pm} r_1^m r_3^{n+9} \right. \\ \left. + \sum_{\mathbf{a}=\pm} \left\{ A_{m,n}^{\pm\pm}(\mathbf{a}) r_1^m r_3^{n-p_3+\frac{13}{2}+i\mathbf{a}\mu_\beta} + B_{m,n}^{\pm\pm}(\mathbf{a}) r_1^{n-p_1-\frac{5}{2}+i\mathbf{a}\mu_\alpha} r_3^{m+9} \right\} \right], \quad (3.23) $$

which comprises three terms, each exhibiting distinct analytic structures about $r_1$ and $r_3$ due to their different origins. The C-term is analytic in both the $r_1$ and $r_3$ directions, indicating that neither of the massive propagators excites CC signals, and thus, it serves purely as the background. On the other hand, the A and B-terms feature oscillatory behavior in one direction, suggesting that one of the massive propagators excites CC signals, making it the signal-background mixing term. By substituting the ansatz into (3.20) and (3.21), we obtain the recurrence relations for the different terms which are summarized at the end of this subsection.

With the explicit forms of the coefficients and after performing one of the infinite summations in (3.23), we finally obtain the particular solution as

$$ \begin{aligned} \hat{\mathcal{P}}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3} = & \frac{u^{p_1+4}}{v^{p_1+5}} \frac{e^{\mp i \frac{\pi}{2} p_{123}}}{2^{10+p_{123}} [\mu_{\alpha}^2+(p_1+5/2)^2]^2} \Gamma\left[\frac{5}{2}+p_3-i\mu_{\beta}, \frac{5}{2}+p_3+i\mu_{\beta}\right] \\ & \times \sum_{m=0}^{\infty} \frac{(3+p_1)_m r_1^m}{\left(\frac{7}{2}+p_1-i\mu_{\alpha}\right)_m \left(\frac{7}{2}+p_1+i\mu_{\alpha}\right)_m} \left\{ r_3^{9+p_{123}} {}_3F_2\left[\frac{1,3+p_3,9+m+p_{123}}{2} \middle| r_3 \right] \right. \\ & \left. \pm \sum_{\mathbf{a}=\pm} (\coth(\pi a \mu_{\beta}) \mp 1) r_3^{\frac{13}{2}+p_{12}-i a \mu_{\beta}} {}_2F_1\left[\frac{1}{2}-i a \mu_{\beta}, \frac{13}{2}+m+p_{12}-i a \mu_{\beta} \middle| r_3 \right] \right\} \\ & + \left( \begin{matrix} \mu_{\alpha} \leftrightarrow \mu_{\beta} \\ p_1 \leftrightarrow p_3 \\ u \leftrightarrow v \end{matrix} \right). \end{aligned} \quad (3.24) $$
The summary of recurrence relations and solutions of different coefficients.

* The coefficient related to the background term $C_{m,n}^{\pm\pm}$:

$$C_{m+1,n}^{\pm\pm} = \frac{(m+3+p_1)(m+n+p_{123}+9)}{\mu_{\alpha}^2+(p_1+m+7/2)^2} C_{m,n}^{\pm\pm}, \quad (3.25)$$

$$C_{0,n}^{\pm\pm} = \frac{\mp i 2 \sqrt{\pi} (3+p_3)_n \Gamma(p_{123}+9+n)}{[\mu_{\alpha}^2+(p_1+5/2)^2][\mu_{\beta}^2+(p_3+5/2)^2] (\frac{7}{2}+p_3-i\mu_{\beta})_n (\frac{7}{2}+p_3+i\mu_{\beta})_n}, \quad (3.26)$$

$$C_{m,n+1}^{\pm\pm} = \frac{(n+3+p_3)(m+n+p_{123}+9)}{\mu_{\beta}^2+(p_3+n+7/2)^2} C_{m,n}^{\pm\pm}, \quad (3.27)$$

$$C_{m,0}^{\pm\pm} = C_{0,n}^{\pm\pm} \quad \text{with} \quad (p_1 \leftrightarrow p_3, \alpha \leftrightarrow \beta, n \rightarrow m). \quad (3.28)$$

* The coefficient related to the background-signal mixing term $A_{m,n}^{\pm\pm}$:

$$A_{m+1,n}^{\pm\pm(a)} = \frac{(m+3+p_1)(m+n+p_{12}+\frac{13}{2}+i\alpha\mu_{\beta})}{\mu_{\alpha}^2+(p_1+m+7/2)^2} A_{m,n}^{\pm\pm(a)}, \quad (3.29)$$

$$A_{m,n+1}^{\pm\pm(a)} = \frac{(n+\frac{1}{2}+i\alpha\mu_{\beta})(m+n+p_{12}+\frac{13}{2}+i\alpha\mu_{\beta})}{(n+1)(n+1+2i\alpha\mu_{\beta})} A_{m,n}^{\pm\pm(a)}, \quad (3.30)$$

$$A_{0,n}^{\pm\pm(a)} = \frac{i e^{\pm\pi a \mu_{\beta}} \sqrt{\pi} \operatorname{csch}(\pi a \mu_{\beta})}{\mu_{\alpha}^2 + (p_1 + 5/2)^2} \times \Gamma\left[ \frac{1}{2} + n + i\alpha\mu_{\beta}, \frac{13}{2} + n + p_{12} + i\alpha\mu_{\beta}, \frac{5}{2} + p_3 - i\mu_{\beta}, \frac{5}{2} + p_3 + i\mu_{\beta} \middle| 1+n, 3+p_3, 1+n+2i\alpha\mu_{\beta} \right]. \quad (3.31)$$

As for the other coefficient $B_{m,n}^{\pm\pm(a)}$, it satisfies the same relation as $A_{m,n}^{\pm\pm(a)}$ with the replacement $p_1 \leftrightarrow p_3$ and $\alpha \leftrightarrow \beta$.

* Finally, solutions of those recurrence relation are given by

$$C_{m,n}^{\pm\pm} = \frac{\mp i 2 \sqrt{\pi} \Gamma(p_{123}+9+m+n)}{[\mu_{\alpha}^2+(p_1+5/2)^2][\mu_{\beta}^2+(p_3+5/2)^2]} \times \frac{(3+p_1)_m (3+p_3)_n}{(\frac{7}{2}+p_1-i\mu_{\alpha})_m (\frac{7}{2}+p_1+i\mu_{\alpha})_m (\frac{7}{2}+p_3-i\mu_{\beta})_n (\frac{7}{2}+p_3+i\mu_{\beta})_n}, \quad (3.32)$$

$$A_{m,n}^{\pm\pm(a)} = \frac{i\sqrt{\pi} e^{\pm\pi a \mu_{\beta}} \operatorname{csch}(\pi a \mu_{\beta})}{\mu_{\alpha}^2+(p_1+5/2)^2} \frac{(3+p_1)_m}{(\frac{7}{2}+p_1-i\mu_{\alpha})_m (\frac{7}{2}+p_1+i\mu_{\alpha})_m} \times \Gamma\left[\frac{5}{2}+p_3-i\mu_{\beta}, \frac{5}{2}+p_3+i\mu_{\beta}, m+n+p_{12}+\frac{13}{2}+i\alpha\mu_{\beta}, \frac{1}{2}+n+i\alpha\mu_{\beta} \middle| 1+n, 3+p_3, 1+n+2i\alpha\mu_{\beta}\right], \quad (3.33)$$

$$B_{m,n}^{\pm\pm(a)} = A_{m,n}^{\pm\pm(a)} \quad \text{with} \quad (p_1 \leftrightarrow p_3, \alpha \leftrightarrow \beta). \quad (3.34)$$
## 3.2 Determine the coefficients

Summarizing the results so far and recalling the relation $\mathcal{I} = u^{-4-p_1}v^{-4-p_3}\hat{\mathcal{I}}$, the general solutions for the seed integral $\mathcal{I}$ are given by

$$ \mathcal{I}_{abc,\alpha\beta}^{p_1 p_2 p_3} = \sum_{d,e=\pm} c_{abc,de} u^{-\frac{5}{2}-p_1-id\mu_\alpha} v^{-\frac{5}{2}-p_3-ie\mu_\beta} \mathcal{F}_4 \left[ \frac{4+p_2-i(d\mu_\alpha+e\mu_\beta)}{2} \middle| \frac{5+p_2-i(d\mu_\alpha+e\mu_\beta)}{2} \right] u^2, v^2 + \frac{\hat{\mathcal{P}}_{abc,\alpha\beta}^{p_1 p_2 p_3}}{u^{p_1+4}v^{p_3+4}}, \quad (3.35) $$

where the particular solutions are given by $\hat{\mathcal{P}}_{\pm\mp,\alpha\beta}^{p_1 p_2 p_3} = 0$, (3.16) and (3.24) respectively. $c_{abc,de}$ are undetermined coefficients and should be fixed by boundary conditions. Previous works on single-exchange diagrams typically use the folded limit and factorised limit as boundary conditions [106], but obtaining the formula under such limits can be challenging in our case. Instead, we adopt an alternative method proposed in [135], where the double soft limit $u, v \to 0$ is imposed as the boundary conditions. Under this limit, the seed integral (3.35) becomes

$$ \lim_{u,v \ll 1} \mathcal{I}_{abc,\alpha\beta}^{p_1 p_2 p_3} = \sum_{d,e=\pm} c_{abc,de} u^{-\frac{5}{2}-p_1-id\mu_\alpha} v^{-\frac{5}{2}-p_3-ie\mu_\beta} \Gamma \left[ \frac{4+p_2-i(d\mu_\alpha+e\mu_\beta)}{2} \middle| \frac{5+p_2-i(d\mu_\alpha+e\mu_\beta)}{2} \right]. \quad (3.36) $$

Note that the particular solutions (the second term of (3.35)) are subdominant in this limit. On the other hand, the time integral under such limit can be readily determined by employing the partial Mellin-Barnes representation, as detailed in the appendix C. By comparing the expressions, we can fix the coefficients as

$$ c_{\pm\mp,de} = -e^{\mp i \frac{\pi}{2} (p_{13} - p_2)} \text{csch}(\pi d \mu_{\alpha}) \text{csch}(\pi e \mu_{\beta}) \times \tilde{\Gamma}(p_1, p_2, p_3, \mu_{\alpha}, \mu_{\beta}), \quad (3.37) $$

$$ c_{\pm\mp,de} = \mp i e^{\mp i \frac{\pi}{2} (p_{12} - p_3)} \text{csch}(\pi d \mu_{\alpha}) \text{csch}(\pi e \mu_{\beta}) e^{\mp \pi d \mu_{\alpha}} \times \tilde{\Gamma}(p_1, p_2, p_3, \mu_{\alpha}, \mu_{\beta}), \quad (3.38) $$

$$ c_{\pm\mp,de} = c_{\mp\mp,ed} \quad \text{with} \quad (p_1 \leftrightarrow p_3) \quad \text{and} \quad (\alpha \leftrightarrow \beta), \quad (3.39) $$

$$ c_{\pm\pm,de} = e^{\mp i \frac{\pi}{2} p_{123}} \text{csch}(\pi d \mu_{\alpha}) \text{csch}(\pi e \mu_{\beta}) e^{\mp \pi(d \mu_{\alpha} + e \mu_{\beta})} \times \tilde{\Gamma}(p_1, p_2, p_3, \mu_{\alpha}, \mu_{\beta}), \quad (3.40) $$

with the new defined function,

$$ \tilde{\Gamma}(p_1, p_2, p_3, \mu_\alpha, \mu_\beta) = \frac{\pi^{\frac{1}{2}}}{2^{4+p_{13}-p_2}} \Gamma \left[ \frac{5}{2} + p_1 - i\mu_\alpha, \frac{5}{2} + p_1 + i\mu_\alpha, \frac{5}{2} + p_3 - i\mu_\beta, \frac{5}{2} + p_3 + i\mu_\beta \middle| 3 + p_1, 3 + p_3 \right]. \quad (3.41) $$

## 3.3 Summary of solutions

After lengthy calculations, we have finally arrived at the exact solution for the double massive exchange seed integral (2.10). In this section, we summarize the detailed forms of various types of seed integrals for reference. We also remind readers again of our notation. The numbers $p_i$ characterize the types of vertices in the seed (2.10), and the summation is denoted by $p_{ij} = p_i + p_j$, etc. The indices $\alpha, \beta$ label the massive fields with masses $m_\alpha$ and $m_\beta$. The momentum ratio is defined as $u = k_1/k_{24}$ and $v = k_3/k_{24}$, where $k_{24} = k_2 + k_4$, and we further introduced $r_1 = 2u/(1+u+v) = 2k_1/k_{1234}$ and $r_3 = 2v/(1+u+v) = 2k_3/k_{1234}$.
The functions ${}_3\mathcal{F}_2$ and $\mathcal{F}_4$ represent the dressed versions of the hypergeometric function and the Appell series, defined in (A.4) and (A.9), respectively. The coefficient $\tilde{\Gamma}(p_1, p_2, p_3, \mu_\alpha, \mu_\beta)$ is defined in eq. (3.41).

* The fully factorised seed integral:

$$ \mathcal{I}_{\pm\pm\pm,\alpha\beta}^{p_1 p_2 p_3} = -e^{\mp i \frac{\pi}{2} (p_{13} - p_2)} \tilde{\Gamma}(p_1, p_2, p_3, \mu_\alpha, \mu_\beta) \sum_{a,b=\pm} \operatorname{csch}(\pi a \mu_\alpha) \operatorname{csch}(\pi b \mu_\beta) u^{-\frac{5}{2}-p_1-i a \mu_\alpha} v^{-\frac{5}{2}-p_3-i b \mu_\beta} \\ \times \mathcal{F}_4 \left[ \frac{4+p_2-i(a\mu_\alpha+b\mu_\beta)}{2}, \frac{5+p_2-i(a\mu_\alpha+b\mu_\beta)}{2} \middle| u^2, v^2 \right]. \quad (3.42) $$

* The partially factorised partially nested seed integral:

$$ \begin{aligned} \mathcal{I}_{\pm\pm\pm,\alpha\beta}^{p_1 p_2 p_3} = & \mp i e^{\mp i \frac{\pi}{2} (p_{12} - p_3)} \tilde{\Gamma}(p_1, p_2, p_3, \mu_\alpha, \mu_\beta) \sum_{a,b=\pm} \operatorname{csch}(\pi a \mu_\alpha) \operatorname{csch}(\pi b \mu_\beta) e^{\mp \pi a \mu_\alpha} u^{-\frac{5}{2}-p_1-i a \mu_\alpha} v^{-\frac{5}{2}-p_3-i b \mu_\beta} \\ & \times \mathcal{F}_4 \left[ \frac{4+p_2-i(a\mu_\alpha+b\mu_\beta)}{2}, \frac{5+p_2-i(a\mu_\alpha+b\mu_\beta)}{2} \middle| u^2, v^2 \right] \\ & + \frac{i}{\pi^{1/2}} \frac{r_3^{p_{12}+\frac{13}{2}}}{v^{p_{12}+3}} \frac{e^{\mp i \frac{\pi}{2} (p_{12}-p_3)}}{2^{6+2p_2}} \tilde{\Gamma}(p_1, p_2, p_3, \mu_\alpha, \mu_\beta) \times \sum_{a=\pm} \sum_{m=0}^{\infty} \Gamma \left[ \frac{3+m+p_1}{\frac{7}{2}+m+p_1-i\mu_\alpha}, \frac{7}{2}+m+p_1+i\mu_\alpha \right] \\ & \times \operatorname{csch}(\pi a \mu_\beta) r_1^m r_3^{i a \mu_\beta} \times {}_2\mathcal{F}_1 \left[ \frac{1}{2}+i a \mu_\beta, \frac{13}{2}+m+p_{12}+i a \mu_\beta \middle| r_3 \right]. \end{aligned} \quad (3.43) $$

The results for $\mathcal{I}_{\pm\pm\pm,\alpha\beta}^{p_1 p_2 p_3}$ can be obtained by replacing $u \leftrightarrow v$ ($r_1 \leftrightarrow r_3$), $p_1 \leftrightarrow p_3$, and $\alpha \leftrightarrow \beta$ in $\mathcal{I}_{\mp\mp\pm,\alpha\beta}^{p_1 p_2 p_3}$.

* The fully nested seed integral:

$$ \begin{aligned} \mathcal{I}_{\pm\pm\pm,\alpha\beta}^{p_1 p_2 p_3} = & e^{\mp i \frac{\pi}{2} p_{123}} \tilde{\Gamma}(p_1, p_2, p_3, \mu_\alpha, \mu_\beta) \times \sum_{a,b=\pm} \operatorname{csch}(\pi a \mu_\alpha) \operatorname{csch}(\pi b \mu_\beta) e^{\mp \pi (a \mu_\alpha + b \mu_\beta)} u^{-\frac{5}{2}-p_1-i a \mu_\alpha} v^{-\frac{5}{2}-p_3-i b \mu_\beta} \\ & \times \mathcal{F}_4 \left[ \frac{4+p_2-i(a\mu_\alpha+b\mu_\beta)}{2}, \frac{5+p_2-i(a\mu_\alpha+b\mu_\beta)}{2} \middle| u^2, v^2 \right] \\ & + \left\{ \frac{r_3^{p_{12}+\frac{13}{2}}}{v^{p_{12}+9}} \cdot \frac{e^{\mp i \frac{\pi}{2} p_{123}}}{2^{10+p_{123}} [\mu_\alpha^2+(p_1+5/2)^2]} \Gamma \left[ \frac{5}{2}+p_3-i\mu_\beta, \frac{5}{2}+p_3+i\mu_\beta \middle| 3+p_3 \right] \right. \\ & \times \sum_{m=0}^{\infty} \frac{(3+p_1)_m r_1^m}{(\frac{7}{2}+p_1-i\mu_\alpha)_m (\frac{7}{2}+p_1+i\mu_\alpha)_m} \left( r_3^{p_3+\frac{5}{2}} {}_3\mathcal{F}_2 \left[ \frac{1,3+p_3,9+m+p_{123}}{\frac{7}{2}+p_3-i\mu_\beta, \frac{7}{2}+p_3+i\mu_\beta} \middle| r_3 \right] \right. \\ & \left. \pm \sum_{a=\pm} (\coth(\pi a \mu_\beta)\mp 1) r_3^{-i a \mu_\beta} \times {}_2\mathcal{F}_1 \left[ \frac{1}{2}-i a \mu_\beta, \frac{13}{2}+m+p_{12}-i a \mu_\beta \middle| r_3 \right] \right) + \left( \begin{array}{c} \mu_\alpha \leftrightarrow \mu_\beta \\ p_1 \leftrightarrow p_3 \\ u \leftrightarrow v \end{array} \right) \right\}. \end{aligned} \quad (3.44) $$
# 4 Continuation to all kinematic regions

## 4.1 Continuation and resummation

The homogeneous solutions presented in the last section are expressed using a specific type of two-variable special function, namely the Appell function $F_4$. This function converges only when $u + v < 1$. Meanwhile, some parts of the particular solution exhibit divergence when $u + v > 1$. Therefore, the aforementioned results about the trispectrum are only applicable when the two massive propagators are relatively soft. Adding to this concern, the physical regions for the bispectrum are constrained to $u + v > 1$, aligning with the triangle inequality $k_1 + k_3 > k_2$.<fnref n="6" /> In this section, we aim to extend the analytical results to all kinematic regions, ensuring validity even in the bispectrum limit. Without loss of generality, we fix the power indices as $p_1 = p_2 = p_3 = -2$ in this section. The results for alternative values of $p_i$ can be derived by applying appropriate differential operators to the subsequent results [135]. Also, we will mainly focus on three representative integrals, namely $\mathcal{I}_{+-+}$, $\mathcal{I}_{++-}$ and $\mathcal{I}_{+++}$, since other components can be obtained by taking complex conjugations or permuting the variables.

**Seed integral $\mathcal{I}_{+-+}$**. Let us begin with the most straightforward case, denoted as $\mathcal{I}_{+-+}$, which is completely factorised in the time integrals and satisfies the homogeneous bootstrap equations, (2.45) and (2.46). To speed up the convergence of series, we observed that employing the variables $r_i \equiv 2k_i/k_{1234}$ or $\tilde{r}_i \equiv 2k_i/k_{123}$ (for the case of the bispectrum), defined as substitutes for $(u, v)$, proves to be more efficient. Note that the divergence for $u + v > 1$ mentioned above appears in the region $r_1 + r_3 > 1$ for the new variables.

Two-variable hypergeometric series $F_4$, which are present in the homogeneous solution (3.42), can be transformed into another second-type Appell function $F_2$ with arguments $r_i$ using the formula (A.10). Subsequently, $F_2$ can be expanded as a series of hypergeometric functions through (A.11). This allows us to express the original expression as follows,

$$ \mathcal{I}_{+-+}(r_1, r_3) = \sum_{a,b=\pm} \sum_{m=0}^{\infty} \frac{\pi^2}{2^3} \mathrm{csch}(2\pi a \mu_\alpha) \mathrm{csch}(2\pi b \mu_\beta) (r_1+r_3-2)^3 r_1^{-\frac{1}{2}+m+ia\mu_\alpha} r_3^{-\frac{1}{2}+ib\mu_\beta} \\ \times \Gamma\left[ \frac{1}{2} + m + ia\mu_\alpha \middle| 1+m, 1+m+2ia\mu_\alpha \right] {}_2F_1\left[ \frac{1}{2} + ib\mu_\beta, 2+m+ia\mu_\alpha + ib\mu_\beta \middle| r_3 \right]. \quad (4.1) $$

In the obtained result, there are four series corresponding to indices $a, b = \pm$. One can check numerically that each term diverges when the sum of $r_1$ and $r_3$ exceeds one. However, intriguingly, cancellations occur when all terms are summed together. To better visualize these cancellations, we can regroup the terms in pairs and employ the hypergeometric function

<footnote n="6">⁶A similar divergence in three-point correlators also emerges in the conformal correlation functions in momentum space, specifically in the Triple-K integrals [176].</footnote>
transformation (A.5) to rewrite the expression as

$$ \begin{aligned} \mathcal{I}_{+-+}(r_1, r_3) = & \sum_{m=0}^{\infty} i \pi^2 \operatorname{csch}(2\pi\mu_\alpha) \operatorname{sech}(\pi\mu_\beta) \left(\frac{r_1}{2} + \frac{r_3}{2} - 1\right)^3 \Gamma\left[\frac{1}{2} + m + i\mu_\alpha, 2+m+i\mu_\alpha-i\mu_\beta \middle| 1+m, 1+m+2i\mu_\alpha, \frac{1}{2}+i\mu_\beta\right] \\ & \times r_1^{-\frac{1}{2}+m+i\mu_\alpha} r_3^{-\frac{1}{2}+i\mu_\beta} {}_2F_1\left[\frac{1}{2} + i\mu_\beta, 2+m+i\mu_\alpha+i\mu_\beta \middle| 1-r_3\right] + (\mu_\alpha \rightarrow -\mu_\alpha). \end{aligned} \quad (4.2) $$

The final expression of $\mathcal{I}_{+-+}$ is now well-behaved even in regions where $r_1 + r_3 > 1$, and thus the validity of the analytical expression is extended to all physical kinematic regions. Regarding the bispectrum limit, one can simply replace every $r_i$ above with $\tilde{r}_i$, and the convergence is also guaranteed. Alternatively, we can take a different perspective on the continuation procedure. In the initial steps of solving homogeneous equations, we can directly introduce an ansatz, similar to our approach to finding particular solutions. After the resummation of one variable, the ansatz for the homogeneous solution will take the form outlined above.

**Seed integral $\mathcal{I}_{++-}$.** The case of the partially factorised partially nested seed integral $\mathcal{I}_{++-}$ closely parallels our previous discussion. The particular solution after summing over the series of $r_3^n$ shows great convergence within the three-point region ($k_1 + k_3 > k_2$). Then, the only thing we need to be concerned about is the homogeneous parts $\mathcal{Y}_{++-}$, which can be addressed using the same procedure as before. After changing variables and expanding the results as the hypergeometric series, we noticed that the expressions can be systematically regrouped in pairs, which finally reads as

$$ \begin{aligned} \mathcal{Y}_{++-}(r_1, r_3) = & \sum_{m=0}^{\infty} \left\{ -\frac{\pi^2}{2} [\operatorname{csch}(\pi\mu_\alpha) + \operatorname{sech}(\pi\mu_\alpha)] \operatorname{sech}(\pi\mu_\beta) \Gamma\left[\frac{1}{2} + m + i\mu_\alpha, 2+m+i\mu_\alpha+i\mu_\beta \middle| 1+m, 1+m+2i\mu_\alpha, \frac{1}{2}-i\mu_\beta\right] \right. \\ & \left. \times r_1^{-\frac{1}{2}+m+i\mu_\alpha} r_3^{-\frac{1}{2}-i\mu_\beta} \left(\frac{r_1}{2} + \frac{r_3}{2} - 1\right)^3 {}_2F_1\left[\frac{1}{2}-i\mu_\beta, 2+m+i\mu_\alpha-i\mu_\beta \middle| 1-r_3\right] + (\mu_\alpha \rightarrow -\mu_\alpha) \right\}. \end{aligned} \quad (4.3) $$

Adding the particular solution (3.16), we obtain the final results valid for all physical kinematic regions,

$$ \begin{aligned} \mathcal{I}_{++-}(r_1, r_3) = & \sum_{m=0}^{\infty} \left\{ -\frac{\pi^2}{2} [\operatorname{csch}(\pi\mu_\alpha) + \operatorname{sech}(\pi\mu_\alpha)] \operatorname{sech}(\pi\mu_\beta) \Gamma\left[\frac{1}{2} + m + i\mu_\alpha, 2+m+i\mu_\alpha+i\mu_\beta \middle| 1+m, 1+m+2i\mu_\alpha, \frac{1}{2}-i\mu_\beta\right] \right. \\ & \left. \times r_1^{-\frac{1}{2}+m+i\mu_\alpha} r_3^{-\frac{1}{2}-i\mu_\beta} \left(\frac{r_1}{2} + \frac{r_3}{2} - 1\right)^3 {}_2F_1\left[\frac{1}{2}-i\mu_\beta, 2+m+i\mu_\alpha-i\mu_\beta \middle| 1-r_3\right] + (\mu_\alpha \rightarrow -\mu_\alpha) \right\} \\ & + \sum_{m=0}^{\infty} \left\{ i \pi^2 \operatorname{csch}(2\pi\mu_\beta) \operatorname{sech}(\pi\mu_\alpha) \Gamma\left[\frac{1}{2} + m \middle| \frac{3}{2} + m - i\mu_\alpha, \frac{3}{2} + m + i\mu_\alpha\right] \right. \\ & \left. \times r_1^m r_3^{-\frac{1}{2}+i\mu_\beta} \left(\frac{r_1}{2} + \frac{r_3}{2} - 1\right)^3 {}_2F_1\left[\frac{1}{2} + i\mu_\beta, \frac{5}{2} + m + i\mu_\beta \middle| r_3\right] + (\mu_\beta \rightarrow -\mu_\beta) \right\}, \end{aligned} \quad (4.4) $$
where the first and second summations represent homogeneous and particular solution respectively. Before moving to the discussion about the most challenging case, let us first comment on the order of summation. Initially, the particular solution involves a double summation series over $r_1^m r_3^n$, and we are able to sum up one layer of an infinite series. As we chose earlier, the initial step involves summation over $r_3^n$, resulting in a convergent series. However, an alternative approach involves initially summing over the series related to $r_1$. In this case, the particular solution then becomes

$$ \mathcal{P}_{++-}(r_1, r_3) = \sum_{n=0}^{\infty} i \pi^2 \operatorname{csch}(2\pi\mu_\alpha) \operatorname{sech}(\pi\mu_\alpha) \left(\frac{r_1}{2} + \frac{r_3}{2} - 1\right)^3 \times \Gamma\left[\frac{1}{2} + n + i\mu_\beta \middle| 1 + n, 1 + n + 2i\mu_\beta\right] \\ \times r_3^{-\frac{1}{2}+n+i\mu_\beta} {}_3F_2\left[\begin{array}{c} 1, 1, \frac{5}{2} + n + i\mu_\beta \\ \frac{3}{2} - i\mu_\alpha, \frac{3}{2} + i\mu_\alpha \end{array} \middle| r_1\right] + (\mu_\beta \rightarrow -\mu_\beta), \quad (4.5) $$

but one can verify numerically that this series is divergent for regions where $r_1+r_3 > 1$ (more precisely, it behaves as an asymptotic series). To ensure the final contribution is convergent, instead of expanding it as a series of hypergeometric functions with the argument $r_3$ as in (4.1), we should expand it with argument dependence on $r_1$, which is

$$ \mathcal{Y}_{++-}(r_1, r_3) = \sum_{a,b=\pm} \sum_{m=0}^{\infty} -i \pi^2 e^{\pi a \mu_\alpha} \operatorname{csch}(2\pi a \mu_\alpha) \operatorname{csch}(2\pi b \mu_\beta) \left(\frac{r_1}{2} + \frac{r_3}{2} - 1\right)^3 r_1^{-\frac{1}{2}+i a \mu_\alpha} r_3^{-\frac{1}{2}+m+i b \mu_\beta} \\ \times \Gamma\left[\begin{array}{c} \frac{1}{2}+m+i b \mu_\beta \\ 1+m, 1+m+2 i b \mu_\beta \end{array}\right] {}_2F_1\left[\begin{array}{c} \frac{1}{2}+i a \mu_\alpha, 2+m+i a \mu_\alpha+i b \mu_\beta \\ 1+2 i a \mu_\alpha \end{array}\middle| r_1\right]. \quad (4.6) $$

Unfortunately, regrouping it in pairs, as in (4.3), to achieve convergence is not possible in this expression since the coefficient of each term does not match. This is expected, as the particular solution (4.5) is superficially divergent and requires another component to cancel it. We have explicitly verified that although both (4.5) and (4.6) exhibit superficial divergences at kinematic regions relevant for three-point correlators, the summation of these two is actually well convergent in all physical regions, yielding the same results as (4.4). In other words, we should treat (4.5) and (4.6) as a whole, which provides a complete and finite result. Let us highlight the difference in these two final expressions: for $\mathcal{I}_{++-}$, the argument inside the hypergeometric function is $r_3$, representing the summation over the $r_3$ direction being performed first. On the other hand, for $\mathcal{P}_{++-}$ and $\mathcal{Y}_{++-}$ in equations (4.5) and (4.6), we summed over $r_1$ first. In summary, to obtain the final convergent result, we observe that the order of summation is crucial, and both the particular and homogeneous solutions should be summed along the same direction to eliminate the superficial divergence. Otherwise, we may end up with an unphysical result.

**Seed integral $\mathcal{I}_{+++}$.** Finally, let's examine the fully nested time integral $\mathcal{I}_{+++}$, which poses the most challenging case. By employing the same procedure of expanding the original homogeneous solution, we reach the following result,

$$ \mathcal{Y}_{+++}(r_1, r_3) = \sum_{a,b=\pm} \sum_{m=0}^{\infty} \pi^2 e^{\pi(a\mu_\alpha + b\mu_\beta)} \operatorname{csch}(2\pi a\mu_\alpha) \operatorname{csch}(2\pi b\mu_\beta) \left(\frac{r_1}{2} + \frac{r_3}{2} - 1\right)^3 r_1^{-\frac{1}{2}+m+i a \mu_\alpha} r_3^{-\frac{1}{2}+i b \mu_\beta} \\ \times \Gamma\left[\begin{array}{c} \frac{1}{2}+m+i a \mu_\alpha \\ 1+m, 1+m+2 i a \mu_\alpha \end{array}\right] {}_2F_1\left[\begin{array}{c} \frac{1}{2}+i b \mu_\beta, 2+m+i a \mu_\alpha+i b \mu_\beta \\ 1+2 i b \mu_\beta \end{array}\middle| r_3\right]. \quad (4.7) $$
This homogeneous component remains superficially divergent at kinematic regions relevant for three-point correlators and cannot be regrouped into pairs to cancel out the divergence. Similar to the preceding discussion, we therefore need to incorporate contributions from both the homogeneous and particular solutions to obtain finite and convergent results. By incorporating the particular solution, the final expression takes the form,

$$
\begin{equation}
\begin{aligned}
\mathcal{I}_{+++}(r_1, r_3) = & \sum_{\mathbf{a},\mathbf{b}=\pm} \sum_{m=0}^{\infty} \pi^2 e^{\pi(\mathbf{a}\mu_\alpha + \mathbf{b}\mu_\beta)} \mathrm{csch}(2\pi\mathbf{a}\mu_\alpha) \mathrm{csch}(2\pi\mathbf{b}\mu_\beta) \left(\frac{r_1}{2} + \frac{r_3}{2} - 1\right)^3 r_1^{-\frac{1}{2}+m+\mathrm{i}\mathbf{a}\mu_\alpha} r_3^{-\frac{1}{2}+\mathrm{i}\mathbf{b}\mu_\beta} \\
& \times \Gamma\left[\frac{1}{2}+m+\mathrm{i}\mathbf{a}\mu_\alpha \quad 1+m, 1+m+2\mathrm{i}\mathbf{a}\mu_\alpha\right] {}_2F_1\left[\frac{1}{2}+\mathrm{i}\mathbf{b}\mu_\beta, 2+m+\mathrm{i}\mathbf{a}\mu_\alpha+\mathrm{i}\mathbf{b}\mu_\beta \middle| r_3\right] \\
& + \left\{ \sum_{m=0}^{\infty} \frac{\pi^2}{2} [\mathrm{csch}(\pi\mu_\alpha) - \mathrm{sech}(\pi\mu_\alpha)] \mathrm{sech}(\pi\mu_\beta) \left(\frac{r_1}{2} + \frac{r_3}{2} - 1\right)^3 r_1^{-\frac{1}{2}+m-\mathrm{i}\mu_\alpha} \right. \\
& \left. \times \Gamma\left[\frac{1}{2}+m-\mathrm{i}\mu_\alpha \quad 1+m, 1+m-2\mathrm{i}\mu_\alpha\right] {}_3F_2\left[1, 1, \frac{5}{2}+m-\mathrm{i}\mu_\alpha \middle| r_3\right] + (\mu_\alpha \rightarrow -\mu_\alpha) \right\} \\
& + \left\{ \sum_{m=0}^{\infty} \frac{\pi^2}{2} [\mathrm{csch}(\pi\mu_\beta) - \mathrm{sech}(\pi\mu_\beta)] \mathrm{sech}(\pi\mu_\alpha) \left(\frac{r_1}{2} + \frac{r_3}{2} - 1\right)^3 r_1^m r_3^{-\frac{1}{2}-\mathrm{i}\mu_\beta} \right. \\
& \left. \times \Gamma\left[\frac{m+1}{2}+m-\mathrm{i}\mu_\alpha, \frac{3}{2}+m+\mathrm{i}\mu_\alpha\right] {}_2F_1\left[\frac{1}{2}-\mathrm{i}\mu_\beta, \frac{5}{2}+m-\mathrm{i}\mu_\beta \middle| r_3\right] + (\mu_\beta \rightarrow -\mu_\beta) \right\} \\
& + \sum_{m=0}^{\infty} \pi^2 \mathrm{sech}(\pi\mu_\beta) \mathrm{sech}(\pi\mu_\alpha) \left(\frac{r_1}{2} + \frac{r_3}{2} - 1\right)^3 r_1^m \\
& \times \Gamma\left[\frac{1+m}{2}+m-\mathrm{i}\mu_\alpha, \frac{3}{2}+m+\mathrm{i}\mu_\alpha\right] {}_3F_2\left[\frac{1}{2}-\mathrm{i}\mu_\beta, \frac{3}{2}+m+\mathrm{i}\mu_\beta \middle| r_3\right],
\end{aligned}
\tag{4.8}
\end{equation}
$$

where the first summation comes from the homogeneous solution, while the second and third summations arise from the particular solution with a partially non-analytical structure (signal-background mixing term). The last summation is the result of the particular solution, which is analytical in both $r_1$ and $r_3$ (purely background term). As previously emphasized, the order of summation is crucial, and in the formula above, each term involves a summation along $r_3$. Alternatively, one can first sum along the $r_1$ direction, yielding similar results. In figure 3, we compare the exact analytical results (4.8) with the corresponding numerical calculations. We have chosen $r_1 = r_3 = r$, while the mass parameters $\mu_\alpha$ and $\mu_\beta$ are different in the left and right panels. We can see that the analytical expression shows excellent agreement with the numerical calculations. Interestingly, unlike the single-exchange diagram, the oscillatory pattern shown in figure 3 displays novel features, akin to the overlay of two oscillation signals. In section 5, as we delve into phenomenology, we will explore these characteristics in detail.

In summary, certain series may initially exhibit superficial divergences, nevertheless, by properly organizing and treating them as a whole, the final expression yields a finite result applicable to all configurations of physical momentum. It is anticipated that certain terms in the final results can be related to each other through some mathematical formulae, potentially leading to the explicit cancellation of divergent components. For example, one potential approach involves employing the expansion series of the hypergeometric function around the argument unity [177], particularly in the regions where $r_1 + r_3 > 1$. However,
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The comparison between the numerical calculation and exact analytical result (4.8). The gray dots and red solid line represent numerical and analytical results, respectively. For the numerical calculation, we directly perform the time integral of the seed integral (2.10). In both figure we set $r_1 = r_3 = r$. The masses are chosen as $\mu_\alpha = \mu_\beta = 2$ in the left panel and $\mu_\alpha = 1.5, \mu_\beta = 2$ in the right panel. For better visualization, we multiplied $\mathcal{I}_{++}$ by a factor $r^{1/2}$.</figcaption></figure>

this transformation introduces significant complexity and further complicates the expression. These transformations and simplifications are beyond the scope of this paper, and we leave a more detailed discussion to future work. As a final note, we stress that although our final analytical results are still expressed as (single-layer) infinite summations, the series expansions are explicitly convergent. For concrete evaluation, one may therefore incorporate any number of terms to reach an arbitrary precision. For example, in figure 3, we included all terms with $m \le 20$, an approximation which perfectly matches the numerical calculation of the integrals, as far as the eye can tell, and the summations could always be more precise by adding more terms.

## 4.2 Consistency checks at different limits

To further verify the correctness of the analytical expressions, in this subsection, we will perform a detailed examination of their behaviors under different limits. For example, with the standard Bunch-Davies vacuum choice, some spurious poles should cancel out, such that the total results should be smooth when $k_1 = k_{234}$ or $k_3 = k_{124}$, equivalently, $r_1$ or $r_3$ approaches the unity, which is similar to the folded singularity cancellation present in the single-exchange diagrams [106, 178, 179].

**Double soft limits.** The first simple and important limit is the case when both $r_1$ and $r_3$ approach zero, which is the regime in which the non-analytic cosmological collider signals will become the dominant components. This double soft limit is not a consistency check but rather an interesting limit for the phenomenology of the trispectrum which will be discussed in section 5. Under such a limit, the seed integrals become,

$$ \lim_{r_1,r_3 \to 0} \mathcal{I}_{+-}(r_1, r_3) = -\frac{\pi^3 \operatorname{csch}(2\pi\mu_\alpha) \operatorname{csch}(2\pi\mu_\beta)}{2^{2i\mu_\alpha-2i\mu_\beta}} \times \Gamma \left[ \begin{array}{c} 2+i\mu_\alpha-i\mu_\beta \\ 1+i\mu_\alpha, 1-i\mu_\beta \end{array} \right] r_1^{-\frac{1}{2}+i\mu_\alpha} r_3^{-\frac{1}{2}-i\mu_\beta} + \text{c.c.} \\ \sim \mathcal{O} \left( e^{-\pi(\mu_\alpha+2\mu_\beta)} \right) r_1^{-\frac{1}{2}+i\mu_\alpha} r_3^{-\frac{1}{2}-i\mu_\beta} + \text{c.c.}, \quad (4.9) $$
$$
\begin{equation}
\begin{split}
\lim_{r_1, r_3 \to 0} \mathcal{I}_{++-}(r_1, r_3) &= -i e^{\pi \mu_\alpha} \frac{\pi^3 \operatorname{csch}(2\pi \mu_\alpha) \operatorname{csch}(2\pi \mu_\beta)}{2^{2i\mu_\alpha - 2i\mu_\beta}} \times \Gamma \left[ \begin{array}{l} 2+i\mu_\alpha - i\mu_\beta \\ 1+i\mu_\alpha, 1-i\mu_\beta \end{array} \right] r_1^{-\frac{1}{2}+i\mu_\alpha} r_3^{-\frac{1}{2}-i\mu_\beta} \right. \\
&\qquad \left. \sim \mathcal{O} \left( e^{-2\pi \mu_\beta} \right) r_1^{-\frac{1}{2}+i\mu_\alpha} r_3^{-\frac{1}{2}-i\mu_\beta} \right.
\end{split}
\tag{4.10}
\end{equation}
$$

Without loss of generality, we have assumed $\mu_\beta \ge \mu_\alpha$ and neglected terms with heavier Boltzmann suppression in the second lines. For the fully nested seed integral $\mathcal{I}_{+++}$, let us investigate separately contributions from both the homogeneous solution and the particular solution which is

$$
\begin{equation}
\begin{split}
\lim_{r_1, r_3 \to 0} \mathcal{Y}_{+++}(r_1, r_3) &= -e^{\pi(\mu_\alpha + \mu_\beta)} \frac{\pi^3 \operatorname{csch}(2\pi\mu_\alpha) \operatorname{csch}(2\pi\mu_\beta)}{2^{2i\mu_\alpha + 2i\mu_\beta}} \times \Gamma \left[ \begin{array}{c} 2+i\mu_\alpha + i\mu_\beta \\ 1+i\mu_\alpha, 1+i\mu_\beta \end{array} \right] r_1^{-\frac{1}{2}+i\mu_\alpha} r_3^{-\frac{1}{2}+i\mu_\beta} \\
&\sim \mathcal{O} \left( e^{-\pi(\mu_\alpha + \mu_\beta)} \right) r_1^{-\frac{1}{2}+i\mu_\alpha} r_3^{-\frac{1}{2}+i\mu_\beta},
\end{split}
\tag{4.11}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
\lim_{r_1, r_3 \to 0} \mathcal{P}_{+++}(r_1, r_3) &= \frac{\pi^{3/2} \operatorname{csch}(\pi \mu_\beta) + \operatorname{sech}(\pi \mu_\beta)}{2^{2i\mu_\beta - 1} (1 + 4\mu_\alpha^2)} \times \Gamma \begin{bmatrix} \frac{5}{2} + i\mu_\beta \\ 1 + i\mu_\beta \end{bmatrix} r_3^{-\frac{1}{2} + i\mu_\beta} + \begin{pmatrix} \mu_\alpha \leftrightarrow \mu_\beta \\ r_1 \leftrightarrow r_3 \end{pmatrix} \\
&\quad - \frac{32}{(1 + 4\mu_\alpha^2)(1 + 4\mu_\beta^2)} \\
&\sim \mathcal{O} \left( e^{-\pi \mu_\beta} \right) r_3^{-\frac{1}{2} + i\mu_\beta} + \mathcal{O} \left( e^{-\pi \mu_\alpha} \right) r_1^{-\frac{1}{2} + i\mu_\alpha} + \text{constant.}
\end{split}
\tag{4.12}
\end{equation}
$$

It is clear that each signal is accompanied by a corresponding Boltzmann factor. In the homogeneous part $\mathcal{Y}_{+++}$, both massive fields contribute to the CC signals, resulting in a total amplitude of approximately $\mathcal{O}(e^{-2\pi\mu})$. Regarding the particular solutions, some components exhibit a non-analytic structure along one direction, leading to a suppression of $\mathcal{O}(e^{-\pi\mu})$, while entirely analytic terms remain unaffected by any suppression.

**Double folded limit.** As we mentioned previously, an essential consistency check is to consider the limit where both $r_i$ approach the unity. With the standard Bunch-Davis initial condition, the final result should be finite under this limit. Any spurious divergence present in each term of the analytical expression should cancel out. This is similar to the single massive exchange case, where the spurious divergence under the folded limit is also eliminated [106]. Let us start from the simplest $\mathcal{I}_{+-+}$ case. After regrouping the result as in (4.2), it is then easy to take such a limit, where

$$
\begin{equation}
\begin{split}
\lim_{r_3 \to 1} \mathcal{I}_{++-}(r_1, r_3) &= \sum_{m=0}^{\infty} i \frac{\pi^2}{2} \frac{\operatorname{csch}(2\pi\mu_\alpha)\operatorname{sech}(\pi\mu_\beta)(r_1-1)^3}{(1+2m+2i\mu_\alpha)(3+2m+2i\mu_\alpha)} \Gamma \left[ \begin{array}{l} 2+m+i\mu_\alpha-i\mu_\beta, 2+m+i\mu_\alpha+i\mu_\beta \\ 1+m, 1+m+2i\mu_\alpha \end{array} \right] \times r_1^{-\frac{1}{2}+m+i\mu_\alpha} \\
&\quad + (\mu_\alpha \to -\mu_\alpha) \\
&= -i\pi^2 \operatorname{csch}(2\pi\mu_\alpha) \operatorname{sech}(\pi\mu_\beta) \left( \frac{r_1}{2} - \frac{1}{2} \right)^3 r_1^{-\frac{1}{2}-i\mu_\alpha} \times {}_3F_2 \left[ \begin{array}{c} \frac{1}{2}-i\mu_\alpha, 2-i\mu_\alpha-i\mu_\beta, 2-i\mu_\alpha+i\mu_\beta \\ \frac{5}{2}-i\mu_\alpha, 1-2i\mu_\alpha \end{array} \middle| r_1 \right] \\
&\quad + (\mu_\alpha \to -\mu_\alpha).
\end{split}
\tag{4.13}
\end{equation}
$$

Clearly, there is no divergence as $r_3$ approaches unity. Regarding the limit as $r_1$ tends to one, one simply needs to sum over the $r_1$ series first in equation (4.1) and then take the
limit, resulting in a similar expression. The double-folded limit $r_1, r_3 \to 1$ can also be easily verified based on our results (4.13), which is

$$ \lim_{r_1, r_3 \to 1} \mathcal{I}_{+-+}(r_1, r_3) = - \left( \frac{k_{24}}{k_1} \right)^3 \pi^2 \mu \operatorname{csch}^2(2\pi\mu) \sinh(\pi\mu). \quad (4.14) $$

This is $\mathcal{O}(e^{-3\pi\mu})$ and here for simplicity, we assume equal masses for both massive exchanges ($\mu_\alpha = \mu_\beta = \mu$). The prefactor originates from $(r_1 + r_3 - 2)^3$ in the initial expression (4.2), and we have used the relation $r_i = 2k_i/k_{1234}$. As for the seed integral $\mathcal{I}_{++}$, through the same procedural steps, we observe the cancellation of spurious poles upon summing all contributions. Consequently, the finite terms under the double folded limit $r_1, r_3 \to 1$ become

$$ \lim_{r_1, r_3 \to 1} \mathcal{I}_{++-}(r_1, r_3) = - \left( \frac{k_{24}}{k_1} \right)^3 \left[ \frac{\pi^2}{32} (4\mu^2 - 3) \operatorname{sech}^2(\pi\mu) + i \frac{\pi^2 \mu}{2} \operatorname{csch}(2\pi\mu) \right], \quad (4.15) $$

which is $\mathcal{O}(e^{-2\pi\mu})$. The other partially factorised partially nested integral $\mathcal{I}_{-++}$ behaves in the same way. For the fully nested part $\mathcal{I}_{+++}$ (4.8), it can be verified that the poles under the limit $r_3 \to 1$ cancel each other out order by order in the summation. However, expressing the general form for the finite terms is challenging, making it difficult to derive a closed formula under the double folded limit. Instead, as an illustrative example and for the purpose of subsequent application, we can consider the limit relevant for three-point correlators. Then, under the squeezed limit, $\tilde{r}_3 \to 1$ and $\tilde{r}_1 \to 0$, therefore we only need to keep the first order in the summation (4.8) as we will demonstrate below.

**Total energy limit.** Interesting information lies in the limit where the total energy $k_T \equiv k_{1234} \to 0$, as the coefficients of the leading $k_T$ singularities correspond to those of flat spacetime scattering amplitudes [103, 162]. For the diagrams we investigate in this work, the total energy pole is indistinguishable from the partial energy poles due to the linear-mixing between the massive fields $\sigma_\alpha$ and the inflaton fluctuation $\phi$ [150], although this degeneracy can be broken by introducing different sound speeds [131]. The analysis of the leading singularity behavior of the double exchange diagram under a vanishing total energy limit is not as straightforward as in the single-exchange case. Under $k_T \to 0$ limit, the variables we used, $r_1, r_3$, approach infinity. In this limit, the leading singularity arises from the fully nested seed integral $\mathcal{I}_{+++}$ (more precisely, the last term in (4.8), since it is fully connected in the bulk time integral.) Since we have already summed the solution series in one direction and packaged it as some hypergeometric functions with arguments either $r_1$ or $r_3$, then we can first use the expansion of the hypergeometric function as $r_3 \to \infty$, where the leading singularity behavior exhibits some logarithm dependence:

$$ \lim_{k_T \to 0} \mathcal{I}_{+++} = \lim_{k_T \to 0} -\frac{\pi}{8\cosh(\pi\mu)} \frac{r_{13}^3 \log(r_3)}{r_3} \sum_{m=0}^{\infty} \Gamma \left[ \frac{3}{2} + m - i\mu, \frac{3}{2} + m + i\mu \right] r_1^m. \quad (4.16) $$

Here, we have assumed for simplicity equal masses for both massive exchanges. Now, another series with summation variable $m$ can be performed explicitly,

$$ \lim_{k_T \to 0} \mathcal{I}_{+++} = \lim_{k_T \to 0} -\frac{1}{(\mu^2 + 1/4)} \frac{r_{13}^3 \log(r_3)}{8r_3} \times {}_3F_2 \left[ \frac{1, 1, 2}{3} - i\mu, \frac{3}{2} + i\mu \middle| r_1 \right], \quad (4.17) $$
and by expanding around infinity again, we can determine the leading singularity as follows:

$$ \lim_{k_T \to 0} \mathcal{I}_{+++} = \frac{r_{13}^3}{8 r_1 r_3} \log(r_1) \log(r_3) \sim \frac{k_{13}^3}{4k_1 k_3 k_T} \log^2(k_T). \quad (4.18) $$

Thus, we find that the coefficient is independent of the masses. This was expected since the $k_T$ poles arise from the past infinity $\tau \to -\infty$ where all mass terms become negligible. This form of the leading singularity (4.18) can also be found using the early time expansion of the mode functions without the knowledge of the full solution. For example, the mixing propagator between $\sigma$ and $\phi$, see [131], under the limit is:

$$ \lim_{\tau \to -\infty} \mathcal{K}(k, \tau) \equiv \lim_{\tau \to -\infty} \int_{-\infty}^{0} d\tau' (\tau')^{-2} e^{ik\tau'} D_{++}(k; \tau, \tau') \sim -\frac{\tau e^{ik\tau}}{2k} \log(-k\tau), \quad (4.19) $$

so the fully nested integral under such limit becomes

$$ \begin{aligned} \mathcal{I}_{+++} = & -i k_{24}^3 \int_{-\infty}^{0} d\tau \frac{e^{ik_{24}\tau}}{\tau^2} \mathcal{K}(k_1, \tau) \mathcal{K}(k_3, \tau) \sim -i \frac{k_{24}^3}{4k_1 k_3} \int_{-\infty}^{0} d\tau e^{ik_T \tau} \log(-k_1 \tau) \log(-k_3 \tau) \\ & \sim -\frac{k_{24}^3}{4k_1 k_3} \frac{\log^2(k_T)}{k_T}, \end{aligned} \quad (4.20) $$

and noting that $k_{24} = -k_{13}$ under $k_T \to 0$, it exactly matches (4.18).

**Squeezed limit.** Another intriguing limit that will play significant roles in the subsequent discussion is the so-called squeezed limit where $\tilde{r}_3(r_3) \to 1$ and $\tilde{r}_1(r_1) \to 0$. For the fully factorised seed integral, the result under $r_3 \to 1$ is already obtained in (4.13) and we only need to take the series expansion around $r_1 = 0$ which gives

$$ \lim_{\substack{r_3 \to 1 \\ r_1 \to 0}} \mathcal{I}_{+-+}(r_1, r_3) = \frac{\pi^2 \operatorname{sech}^2(\pi\mu) \operatorname{csch}(\pi\mu)}{4(3i - 2\mu)} \times r_1^{-\frac{1}{2}+i\mu} + (\mu \to -\mu), \quad (4.21) $$

with an amplitude $\mathcal{O}(e^{-3\pi\mu})$. For two partially factorised partially nested integrals, following the same way, the results are

$$ \lim_{\substack{r_3 \to 1 \\ r_1 \to 0}} \mathcal{I}_{++-}(r_1, r_3) = \frac{\pi^2 [1 - \tanh(\pi\mu)][1 + \coth(2\pi\mu)]}{2(3 + 2i\mu)} \times r_1^{-\frac{1}{2} + i\mu} + (\mu \to -\mu), \quad (4.22) $$

$$ \lim_{\substack{r_3 \to 1 \\ r_1 \to 0}} \mathcal{I}_{--+}(r_1, r_3) = \left( \frac{\pi^2 \operatorname{sech}^2(\pi\mu) [1 + \coth(\pi\mu)]}{4(2i\mu - 3)} + \frac{\pi^{1/2}}{2^{4-2i\mu}} \frac{2\mu + 5i}{2\mu + 3i} \Gamma[\frac{1}{2} - i\mu, i\mu] \operatorname{sech}(\pi\mu) \right) \times r_1^{-\frac{1}{2} - i\mu} \\ + (\mu \to -\mu), \quad (4.23) $$

where both of them have an amplitude $\mathcal{O}(e^{-2\pi\mu})$. For the fully nested seed integral $\mathcal{I}_{+++}$, as mentioned previously, under $r_1 \to 0$, we will only retain the first order in the series summation. It is straightforward to show that the poles that would diverge as $r_3 \to 1$ are canceled, yielding the final finite contributions

$$ \lim_{\substack{r_3 \to 1 \\ r_1 \to 0}} \mathcal{I}_{+++}(r_1, r_3) = \left( \frac{\pi^2 \operatorname{csch}(2\pi\mu) \operatorname{sech}(\pi\mu)}{2(\mu - 3i)} + \frac{\pi^{1/2}}{2^{4+2i\mu}} \frac{2i\mu + 5}{2\mu - 3i} \Gamma[\frac{1}{2} + i\mu, -i\mu] (1 + \tanh(\pi\mu)) \right) r_1^{-\frac{1}{2} + i\mu} \\ + (\mu \to -\mu), \quad (4.24) $$
where the first term arises from the homogeneous solution and is suppressed by $O(e^{-3\pi\mu})$. In contrast, the second term comes from one of the particular solutions, wherein only one massive propagator contributes to the signal and the other massive line lies within the analytic region and contributes to the background. Hence, the overall magnitude of the second term is only suppressed by $O(e^{-\pi\mu})$.

For the bispectrum in the squeezed limit, the dominant signal arises from the CC signals, with the leading contribution experiencing the least Boltzmann suppression from $\mathcal{I}_{+++}$. Substituting this expression into (2.15), we obtain an approximate result for the three-point correlator in the squeezed limit which is useful for later discussion

$$ \lim_{k_1 \to 0} \langle \varphi_{\mathbf{k}_1} \varphi_{\mathbf{k}_2} \varphi_{\mathbf{k}_3} \rangle' = - \frac{\rho^2 \lambda H}{(k_1 k_2 k_3)^2} \mathrm{Re} \left\{ \left[ \frac{\pi^{1/2}}{2^{4+2i\mu}} \frac{2i\mu+5}{2\mu-3i} \Gamma\left[\frac{1}{2}+i\mu, -i\mu\right] (1+\tanh(\pi\mu)) + \mathcal{O}\left(e^{-2\pi\mu}\right) \right] \left(\frac{k_1}{k_3}\right)^{\frac{1}{2}+i\mu} + \mathcal{O}\left(\frac{k_1}{k_3}\right) \right\}, \quad (4.25) $$

where $\rho_\alpha$ and $\lambda_{\alpha\beta}$ are reduced to $\rho$ and $\lambda$ for the equal mass case, and the $O(k_1/k_3)$ term arises from the leading analytical background contribution. Since $\mathcal{I}_{+++}$ dominates in the squeezed limit, and because we will also discuss the CC signals with multiple isocurvature species in the next section, for later convenience, here we also show the expression when two propagators have different masses explicitly:

$$ \lim_{\substack{r_3 \to 1 \\ r_1 \to 0}} \mathcal{I}_{+++}(r_1, r_3) = -i \frac{\pi^{3/2}}{2^{4+2i\mu_\alpha}} \mathrm{sech}(\pi\mu_\beta) [1+\tanh(\pi\mu_\alpha)] \times \Gamma\left[-1-i\mu_\alpha+i\mu_\beta, -1-i\mu_\alpha-i\mu_\beta\right] \\ \times r_1^{-\frac{1}{2}+i\mu_\alpha} {}_3F_2\left[-\frac{3}{2}-i\mu_\alpha, -1-i\mu_\alpha-i\mu_\beta, -1-i\mu_\alpha+i\mu_\beta; 1\right] + (\mu_\alpha \to -\mu_\alpha), \quad (4.26) $$

and the three-point function in the squeezed limit becomes:

$$ \lim_{k_1 \to 0} \langle \varphi_{\mathbf{k}_1} \varphi_{\mathbf{k}_2} \varphi_{\mathbf{k}_3} \rangle' = \sum_{\alpha, \beta} \frac{\rho_\alpha \rho_\beta \lambda_{\alpha\beta} H}{(k_1 k_2 k_3)^2} \mathrm{Re} \left\{ \left[ i \frac{\pi^{3/2}}{2^{4+2i\mu_\alpha}} \mathrm{sech}(\pi\mu_\beta) [1+\tanh(\pi\mu_\alpha)] \times \Gamma\left[-1-i\mu_\alpha+i\mu_\beta, -1-i\mu_\alpha-i\mu_\beta\right] \right. \right. \\ \left. \left. \times {}_3F_2\left[-\frac{3}{2}-i\mu_\alpha, -1-i\mu_\alpha-i\mu_\beta, -1-i\mu_\alpha+i\mu_\beta; 1\right] + \mathcal{O}(e^{-2\pi\mu_\alpha}, e^{-2\pi\mu_\beta}) \right] \left(\frac{k_1}{k_3}\right)^{\frac{1}{2}+i\mu_\alpha} + \mathcal{O}\left(\frac{k_1}{k_3}\right) \right\}. \quad (4.27) $$

## 5 Double massive exchange phenomenology

In this section, we discuss the physical effects that can be read off our analytical solution, eqs. (4.2), (4.4) and (4.8). We also compare it to the predictions from different techniques, namely the effective single-field description when the massive fields are integrated out and an exact numerical evolution that does not rely on a perturbative scheme in terms of quadratic mixings. In the subsequent sections, we explore different regimes of interest. Before that, we start by presenting the motivation for the setup and these complementary techniques.

**Field content, interactions and relations to multifield inflation.** For almost all concrete applications of this section, we focus on a simple field content made of the inflaton's massless fluctuations $\varphi$ and a single extra massive scalar field $\sigma$ (see an exception in section 5.4).
These are already mixing at the level of the quadratic interactions, as encoded by eq. (2.11) with $\rho_{\alpha} = \rho \times \delta_{\alpha 1}$ and $\sigma_1 = \sigma$. Additional cubic interactions leading to double-exchange diagrams that we consider are described by eq. (2.12) with $\lambda_{\alpha\beta} = \lambda \times \delta_{\alpha 1}\delta_{\beta 1}$ and are leading to a three-point correlator of $\varphi$ given by the double-exchange seed integral in the limit $k_4 \to 0$ and with $(p_1, p_2, p_3) = (-2, -2, -2)$, see eq. (2.15). Considering these particular interactions is motivated by different though complementary approaches to multifield inflation, as well as by the potential of detectability of this channel.

First, popular multifield realizations of the inflationary paradigm, called general non-linear sigma models and encoding both potential and kinetic interactions, generically predict their existences. In this two-field context, the quadratic mixing $\rho$ can be related to the dimensionless rate of turn of the background trajectory in field space, $\eta_{\perp}$, as $\rho = 2H\eta_{\perp}$, while the cubic interaction $\lambda$ is given as a combination of the turn rate and the scalar curvature of the field space of mass dimension $-2$, $R_{\text{fs}}$, as $\lambda = (H^2\eta_{\perp}^2 - \epsilon H^2 M_{\text{Pl}}^2 R_{\text{fs}})/f_{\pi}^2$ [174]. In this expression, $f_{\pi}^2$ represents the normalisation of the curvature fluctuation with respect to the massless inflaton's fluctuations $\varphi = -(f_{\pi}^2/H) \times \zeta$. Therefore, in the weak quadratic mixing regime where the power spectrum of $\varphi$ is negligibly affected by its interactions with the massive fields, one can relate the typical size of the curvature fluctuation $\Delta_{\zeta}$ to the one of a massless uncoupled field in the Sitter $\Delta_{\varphi} = H/(2\pi)$ as

$$ \Delta_{\zeta} \equiv A_s^{1/2} = \frac{H}{f_{\pi}^2} \Delta_{\varphi} = \frac{H^2}{2\pi f_{\pi}^2}, \qquad (5.1) $$

with $A_s = 2.1 \times 10^{-9}$ [180] at CMB scales. This equality enables one to express $f_{\pi}^2$ as a large number in units of $H^2$. Actually, for more than one additional massive field $\sigma_{\alpha}$, the form of the interactions depends on the choice of the basis chosen for these fluctuations, leading to a subtle interplay between flavor and mass bases eigenstates as described in ref. [77]. In the latter reference, the form of the interactions in the mass basis is specified, from the knowledge of the ones in the flavor basis in the context of non-linear sigma models with any number of fields, first derived in ref. [175]. We will explore this more generic situation in section 5.4.

Second, generic arguments from the effective field theory of multifield fluctuations tell us that we have to entertain the possibility of these interactions. In this language, the adiabatic curvature fluctuation $\zeta$ is related by a gauge transformation to the Goldstone boson of broken time diffeomorphisms, $\pi$ [181, 182]. At the linear level and when defining the Goldstone boson in the flat gauge, this gauge transformation reads $\zeta = -H\pi + \dots$, where dots stand for order two and more terms in fluctuations. One can then couple it to either massless [183] or massive [3] additional scalar fields $\sigma^{\alpha}$. In the unitary gauge, the adiabatic degree of freedom is hidden in the spacetime metric fluctuation $\delta g^{00}$, that can couple to these additional fields via couplings of the form $b_{\alpha}^{(n)}(\delta g^{00})^n \sigma^{\alpha}$, $c_{\alpha\beta}^{(n)}(\delta g^{00})^n \sigma^{\alpha}\sigma^{\beta}$, etc. Those couplings generically lead to interactions of the forms that we considered above (amongst other ones) [171], after the Stückelberg trick reintroducing the Goldstone boson via, e.g., $\delta g^{00} \to -2\dot{\pi} + \partial_{\mu}\pi\partial^{\mu}\pi$.

Third, the focus of this work is on double-exchange diagrams, where two massive field propagators are exchanged at tree level via the existence of a 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, 90, 170], 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 (\phi')^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 \phi' \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: } |\lambda| \lesssim 1. \qquad (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<fnref n="7" />

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

<footnote n="7">⁷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.</footnote>
Given the observed amplitude of curvature fluctuations at CMB scales, the bound from naturalness is $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 $\phi$, e.g. in eq. (4.25), 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 $\phi$, in the weak mixing regime where $\Delta_{\phi} = H/(2\pi)$, as follows:

$$ \langle \phi_{\mathbf{k}_1} \phi_{\mathbf{k}_2} \phi_{\mathbf{k}_3} \rangle' = - \frac{H^5}{f_{\pi}^2 (k_1 k_2 k_3)^2} S(k_1, k_2, k_3). \qquad (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_{\text{NL}}$'s parameters:

$$ S \sim f_{\text{NL}} \sim (\Delta_{\zeta})^{-1} \times (\rho/H)^2 \times \lambda, \qquad (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 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 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–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:

$$ \mathcal{L}_{\text{EFT}} = \frac{a^3 \epsilon M_{\text{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} \qquad (5.6) $$

$$ \frac{1}{c_s^2} - 1 = \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}. \qquad (5.7) $$

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

$$ S_{\text{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}, \qquad (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_{\text{eq}}(k_1, k_2, k_3) = \left(\frac{k_1}{k_2} + 5 \text{ perms}\right) - \left(\frac{k_3^2}{k_1 k_2} + 2 \text{ perms}\right) - 2. \qquad (5.9)$$

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

$$\cos(S_1, S_2) = \frac{S_1 \cdot S_2}{\sqrt{(S_1 \cdot S_1)(S_2 \cdot S_2)}} \quad \text{with} \qquad (5.10)$$

$$S_1 \cdot S_2 = \int_{1/2}^{1} dx_2 \int_{1-x_2}^{x_2} dx_3 S_1(k, x_2 \times k, x_3 \times k) S_2(k, x_2 \times k, x_3 \times k), \qquad (5.11)$$

one finds $\cos(S_{\text{EFT}}, S_{\text{eq}}) \approx 0.988$. We thus safely dub its amplitude at equilateral configurations the parameter $f_{\text{NL}}^{\text{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_{\text{NL}}^{\text{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}. \qquad (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^2c_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 figure 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_{\text{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}, \qquad (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 figure 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_{\text{NL}}^{\text{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.

## 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
<figure><img src="image_6.png" /><figcaption>**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 as 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_{\text{NL}}^{\text{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$.</figcaption></figure>

in the previous section. Here we instead investigate the shape dependence on ($x_2, x_3$) of the bispectrum $S(k, x_2k, x_3k)$, 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 figure 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 figure 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 \ge 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].

## 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.25). Translated in terms of the primordial curvature fluctuation bispectrum shape function, it reads:

$$ S_{\text{CC,LO}}^{\text{DE}} = \left(\frac{\rho}{H}\right)^2 \times \lambda(2\pi\Delta_{\zeta})^{-1} \times \text{Re} \left\{ \left[ \frac{\pi^{1/2}}{2^{4+2i\mu}} \frac{2i\mu+5}{2\mu-3i} \Gamma\left[\frac{1}{2}+i\mu, -i\mu\right] (1+\tanh(\pi\mu)) + \mathcal{O}(e^{-2\pi\mu}) \right] \kappa^{\frac{1}{2}+i\mu} \right\}, \quad (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_{\text{CC,LO}}^{\text{DE}} \sim \mathcal{O}(\kappa^{1/2}\mu^{-1/2}e^{-\pi\mu})$. 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_{\text{BG}}^{\text{DE}} \sim \mathcal{O}(\mu^{-4})$ 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.2)-(4.4)-(4.8), 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
<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>
<table><thead><tr><th>Mass value</th><th><i>m</i> = <i>H</i></th><th><i>m</i> = 3<i>H</i>/2</th><th><i>m</i> = 2<i>H</i></th><th><i>m</i> = 5<i>H</i>/2</th><th><i>m</i> = 3<i>H</i></th><th><i>m</i> = 7<i>H</i>/2</th></tr></thead><tbody><tr><td>Cos(<i>S</i>, <i>S</i><sub>eq</sub>)</td><td>0.71</td><td>0.87</td><td><b>0.87</b></td><td><b>0.92</b></td><td><b>0.93</b></td><td><b>0.93</b></td></tr><tr><td>Cos(<i>S</i>, <i>S</i><sub>flat</sub>)</td><td>0.92</td><td>0.84</td><td>0.29</td><td>0.72</td><td>0.69</td><td>0.72</td></tr><tr><td>Cos(<i>S</i>, <i>S</i><sub>orth</sub>)</td><td>-0.66</td><td>-0.48</td><td>0.14</td><td>-0.32</td><td>-0.28</td><td>-0.31</td></tr><tr><td>Cos(<i>S</i>, <i>S</i><sub>loc</sub>)</td><td><b>0.99</b></td><td><b>0.90</b></td><td>0.36</td><td>0.76</td><td>0.75</td><td>0.78</td></tr></tbody></table>

<table><thead><tr><th><i>S</i><sub><i>i</i></sub></th><th><i>S</i><sub>eq</sub></th><th><i>S</i><sub>flat</sub></th><th><i>S</i><sub>orth</sub></th><th><i>S</i><sub>loc</sub></th></tr></thead><tbody><tr><td>Cos(<i>S</i><sub><i>i</i></sub>, <i>S</i><sub>eq</sub>)</td><td><b>1</b></td><td>0.46</td><td>0.01</td><td>0.66</td></tr><tr><td>Cos(<i>S</i><sub><i>i</i></sub>, <i>S</i><sub>flat</sub>)</td><td></td><td><b>1</b></td><td>-0.88</td><td>0.91</td></tr><tr><td>Cos(<i>S</i><sub><i>i</i></sub>, <i>S</i><sub>orth</sub>)</td><td></td><td></td><td><b>1</b></td><td>-0.68</td></tr><tr><td>Cos(<i>S</i><sub><i>i</i></sub>, <i>S</i><sub>loc</sub>)</td><td></td><td></td><td></td><td><b>1</b></td></tr></tbody></table>

**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.

$\kappa^{i\mu}$, the series expansion of the $m$-term in the sum takes the following schematic form: $S_{\text{CC},m}^{\text{DE}} \sim \kappa^{1/2+i\mu+m} (1 + \mathcal{O}(\kappa) + \mathcal{O}(\kappa^2) + \dots)$. 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_{\text{CC}}^{\text{DE}} \sim \kappa^{1/2+i\mu} (1 + \mathcal{O}(\kappa^2) + \dots)$. 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_{\text{BG},m}^{\text{DE}} \sim \kappa^1 (1 + \mathcal{O}(\kappa) + \mathcal{O}(\kappa^2) + \dots)$. 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 figure 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.42)). 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_{\text{CC,LO}}^{\text{SE}}}{S_{\text{CC,LO}}^{\text{SE}}} = \lambda(2\pi\Delta_{\zeta})^{-1} \times \frac{2}{\mu^2} + \mathcal{O}(\mu^{-3}), \qquad (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_{\text{CC,LO}}^{\text{SE}} \sim \mathcal{O}(\kappa^{1/2}\mu^{3/2}e^{-\pi\mu})$ and $S_{\text{BG}}^{\text{SE}} \sim \mathcal{O}(\mu^{-2})$). 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 section 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 figure 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.

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_{\text{eff}}^2 = m^2 + \rho^2$, visible in the frequency of the
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>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$.</figcaption></figure>

oscillations in the cosmological collider signal, now dictated by $\mu_{\text{eff}} = \sqrt{\mu^2 + \rho^2/H^2}$. In figure 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 figure 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.
<figure><img src="image_12.png" /><figcaption>**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_{\text{eff}}$. The bootstrap predictions are obtained by the naive replacement $\mu \to \mu_{\text{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.</figcaption></figure>

## 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:

$$S_{\text{CC,LO}}^{\text{SE}} = \left(\frac{\rho}{H}\right)^2 \times \text{Re} \left[ \kappa^{1/2+i\mu} A_{\text{SE}}(\mu) e^{i\delta_{\text{SE}}(\mu)} \right], \qquad (5.16)$$

$$S_{\text{CC,LO}}^{\text{DE}} = \left(\frac{\rho}{H}\right)^2 \times \lambda(2\pi\Delta_{\zeta})^{-1} \times \text{Re} \left[ \kappa^{1/2+i\mu} A_{\text{DE}}(\mu) e^{i\delta_{\text{DE}}(\mu)} \right]. \qquad (5.17)$$

Here $A_{\text{SE,DE}}$ denote the amplitudes of the CC signals and $\delta_{\text{SE,DE}}$ their phases. In general, the phase of CC signals depends on the mass of exchanged fields as well as the form of the couplings [82]. Nevertheless, for concreteness in this phenomenological discussion, we will fix once and for all the interactions to (2.11) and (2.12) since those are the relevant ones in well-motivated non-linear sigma models of inflation. Then the expression for $A_{\text{DE}}(\mu)e^{i\delta_{\text{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}_{\text{SE}}(\mu) e^{i \delta_{\text{SE}}(\mu)} = -i \frac{\pi^{3/2}}{2^{4+2i\mu}} \frac{\Gamma(7/2 + i\mu)}{(2\mu - i)\Gamma(1 + i\mu)} \tanh(\pi\mu) e^{-\pi\mu} [\coth(\pi\mu) - i \censch(\pi\mu) + 1]^2 . \quad (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'^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 figure 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 figure 9, we instead investigate the dependence of an effective phase $\delta_{\text{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.

**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]:

$$ S_{\text{CC,LO}}^{\text{SE,multi}} = \sum_{\alpha} \left(\frac{\rho_{\alpha}}{H}\right)^2 \text{Re} \left[ \kappa^{1/2 + i\mu_{\alpha}} \mathcal{A}_{\text{SE}}(\mu_{\alpha}) e^{i\delta_{\text{SE}}(\mu_{\alpha})} \right], \quad (5.19) $$
<figure><img src="image_13.png" /><figcaption>**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, $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 \to 1$.</figcaption></figure>

where $A_{SE}$ and $\delta_{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

$$S_{\text{CC,LO}}^{\text{DE,multi}} = \sum_{\alpha,\beta} \frac{\rho_{\alpha}\rho_{\beta}}{H^2} \times \lambda_{\alpha\beta}(2\pi\Delta_{\zeta})^{-1} \times \text{Re} \left[ \kappa^{1/2+i\mu_{\alpha}} A_{\text{DE}}^{\text{multi}}(\mu_{\alpha}, \mu_{\beta}) e^{i\delta_{\text{DE}}^{\text{multi}}(\mu_{\alpha}, \mu_{\beta})} \right], \quad (5.20)$$

where $A_{\text{DE}}^{\text{multi}}$ and $\delta_{\text{DE}}^{\text{multi}}$ are uniquely determined from our analytical formulae (4.27). 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:

$$S_{\text{CC,LO}}^{\text{DE,multi}} = \sum_{\alpha} \left(\frac{\rho_{\alpha}}{H}\right)^2 \text{Re} \left[ \kappa^{1/2+i\mu_{\alpha}} \tilde{A}_{\text{DE}}^{\alpha} e^{i\tilde{\delta}_{\text{DE}}^{\alpha}} \right], \quad (5.21)$$

$$\text{with } \tilde{A}_{\text{DE}}^{\alpha} e^{i\tilde{\delta}_{\text{DE}}^{\alpha}} = \sum_{\beta} \frac{\rho_{\beta}}{\rho_{\alpha}} \times \lambda_{\alpha\beta}(2\pi\Delta_{\zeta})^{-1} \times A_{\text{DE}}^{\text{multi}}(\mu_{\alpha}, \mu_{\beta}) e^{i\delta_{\text{DE}}^{\text{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
<figure><img src="image_14.png" /><figcaption>**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}$.</figcaption></figure>

beyond the scope of this article, we already note an interesting feature. When there is a single isocurvature species, we have seen in figure 9 that the phase is fixed by the frequency of the CC signal once the couplings are given, 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 fixed by its frequency. However, this relation 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, it is very likely that single-exchange diagrams can be ruled out, if the relation 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 figure 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 $\phi$ 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$) = ($\rho \cos(\theta_{12}), \rho \sin(\theta_{12}))$. 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

$$ \langle \varphi_{\mathbf{k}_1} \varphi_{\mathbf{k}_2} \varphi_{\mathbf{k}_3} \varphi_{\mathbf{k}_4} \rangle' = \frac{H^8 (k_{1234}/4)^3}{f_{\pi}^4 (k_1 k_2 k_3 k_4)^3} \mathcal{T}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4). \quad (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:

$$ \mathcal{T} \sim \tau_{\text{NL}} \sim (\rho/H)^2 \times (H^2 \tilde{\lambda}) \times (2\pi\Delta_{\zeta})^{-2}. \quad (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 $\mathbf{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}{\approx} 0.07 \times \rho^2 \tilde{\lambda} \times \frac{f_{\pi}^4}{H^4} \times \frac{1}{\mu^4}. \quad (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 \to 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,

$$ \mathcal{T}_{\text{approx}} = \mathcal{T}_{\text{ss}} + \mathcal{T}_{\text{sb}} + \mathcal{T}_{\text{bb}}. \quad (5.25) $$
<figure><img src="image_15.png" /><figcaption>**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.8) respectively. In the figure we set $k_2 = k_4 = k$ and $k_1 = k_3 = rk$, the masses are chosen as $\mu_\alpha = \mu_\beta = 2$ and the coupling $\bar{\lambda} \cdot \rho^2$ is set to $10^{-2}$. The left region (light red) is dominated by the term $r^{2i\mu}$ while the right region (light blue) is dominated by the oscillation $r^{i\mu}$. For better visualization, we have multiplied the dimensionless trispectrum $T(rk, k, rk, k)$ by a factor $r^{-7/2}$.</figcaption></figure>

<table><caption>**Table 1.** The leading order mass and $r$ dependence of each term in eq. (5.25).</caption><thead><tr><td>different terms</td><td>mass dependence</td><td>r dependence</td></tr></thead><tbody><tr><td>T<sub>ss</sub></td><td>μ<sup>3/2</sup>e<sup>-2πμ</sup></td><td>r<sup>3+2iμ</sup></td></tr><tr><td>T<sub>sb</sub></td><td>μ<sup>3/2</sup>e<sup>-πμ</sup></td><td>r<sup>7/2+iμ</sup></td></tr><tr><td>T<sub>bb</sub></td><td>μ<sup>-4</sup></td><td>r<sup>4</sup></td></tr></tbody></table>

$T_{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, $T_{ss}$ originates from the homogeneous solution $Y_{++}^{2 \ 0 \ -2}$, under the double soft limit, analogous to eq. (4.11) discussed earlier. For $T_{sb}$, the CC signal is generated only from one propagator, while the other contribution arises from the particular solution, similar to eq. (4.12) under the double-soft limit. Additionally, due to permutations, there is another contribution where only one massive propagator becomes soft, approximated by eq. (4.24). The last term $T_{bb}$, represents the leading-order analytical background.

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. $T_{ss}$ excites CC signals in both massive propagators, thus the signal is suppressed by $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+2i\mu}$. In contrast, $T_{sb}$ is suppressed by $O(e^{-\pi\mu})$ and oscillates as $r^{7/2+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 $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 $O(e^{-2\pi\mu}) \times r^{2i\mu+3}$ and $O(e^{-\pi\mu}) \times r^{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^{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. Lastly, as mentioned previously, there exists another type of double-exchange four-point correlator involving two cubic interactions: one $\phi^2\sigma$ and one $\phi\sigma^2$, along with one quadratic mixing. The analytical expression for this trispectrum can be computed using the methods presented in this work, albeit requiring additional significant efforts. Nonetheless, the distinguishing features of the trispectrum discussed in this subsection, which differentiate between single- and double-exchange diagrams, are expected to be quite general and to hold for this other diagram too. Indeed, the double-frequency oscillations arise from the generation of CC signals in both massive propagators when their momenta are relatively soft. This characteristic feature cannot be mimicked by any single-exchange process. Since the excitation of CC signals at each massive propagator is always accompanied by the Boltzmann suppression factor $O(e^{-\pi\mu})$, the leading-order mass suppression of each term depends on the number of massive propagators that excite CC signals. Therefore, the suppression power summarized in table 1 is likely applicable to other double-exchange diagrams, although the exact coefficients would require detailed analytical calculations, which we leave for future studies.

# 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, it is worth delving deeper into understanding the mathematical structures of these boundary bootstrap differential equations involving multiple exchanged fields, such as those explored in [148, 149] for the three-site chain of conformal scalars.<fnref n="8" /> Sixth, although we did mention inflationary flavor oscillations in the context of double-exchange channels and uncovered first interesting unique features, it remains yet to

<footnote n="8">⁸We thank the anonymous referee for pointing this out.</footnote>
conduct a thorough investigation of how to tell different channels apart in practice from (to-be) observations. Seventh, 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. Eighth, 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.

## Acknowledgments

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 under the project code, IBS-R018-D3, 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.

## 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,

$$ \Gamma[z_1, \dots, z_m] \equiv \Gamma(z_1) \dots \Gamma(z_m), \qquad (A.1) $$

$$ \Gamma \begin{bmatrix} z_1, \dots, z_m \\ w_1, \dots, w_n \end{bmatrix} = \frac{\Gamma(z_1) \dots \Gamma(z_m)}{\Gamma(w_1) \dots \Gamma(w_n)}. \qquad (A.2) $$

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

$$ {}_pF_q \begin{bmatrix} a_1, \dots, a_p \\ b_1, \dots, b_q \end{bmatrix} x = \sum_{n=0}^{\infty} \frac{(a_1)_n \cdots (a_p)_n}{(b_1)_n \cdots (b_q)_n} \frac{x^n}{n!}, \qquad (A.3) $$

where $(a)_n \equiv \Gamma[a+n]/\Gamma[a]$ is the Pochhammer symbol.
We also use its “dressed” version, which is defined by

$$ 
\begin{aligned}
p\mathcal{F}_q \left[ \begin{matrix} a_1, \dots, a_p \\ b_1, \dots, b_q \end{matrix} \middle| x \right] &\equiv \Gamma \left[ \begin{matrix} a_1, \dots, a_p \\ b_1, \dots, b_q \end{matrix} \right] p\mathbb{F}_q \left[ \begin{matrix} a_1, \dots, a_p \\ b_1, \dots, b_q \end{matrix} \middle| x \right] \\
&= \sum_{n=0}^{\infty} \Gamma \left[ \begin{matrix} a_1 + n, \dots, a_p + n \\ b_1 + n, \dots, b_q + n \end{matrix} \right] \frac{x^n}{n!}.
\end{aligned} 
\quad (A.4) 
$$

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

$$ 
\begin{aligned}
& \frac{\sin(\pi(c-a-b))}{\pi} {}_2F_1 \left[ \begin{matrix} a, b \\ c \end{matrix} \middle| x \right] \\
&= \frac{1}{\Gamma[c-a, c-b]} {}_2F_1 \left[ \begin{matrix} a, b \\ a+b-c+1 \end{matrix} \middle| 1-x \right] - \frac{(1-x)^{c-a-b}}{\Gamma[c-a, c-b]} {}_2F_1 \left[ \begin{matrix} c-a, c-b \\ c-a-b+1 \end{matrix} \middle| 1-x \right].
\end{aligned} 
\quad (A.5) 
$$

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

$$ F_2 \left[ a \middle| \begin{matrix} b_1, b_2 \\ c_1, c_2 \end{matrix} \middle| x, y \right] = \sum_{m,n=0}^{\infty} \frac{(a)_{m+n} (b_1)_m (b_2)_n}{(c_1)_m (c_2)_n m! n!} x^m y^n, \quad (A.6) $$

$$ F_4 \left[ \begin{matrix} a, b \\ c_1, c_2 \end{matrix} \middle| x, y \right] = \sum_{m,n=0}^{\infty} \frac{(a)_{m+n} (b)_{m+n}}{(c_1)_m (c_2)_n m! n!} x^m y^n, \quad (A.7) $$

and their dressed versions are expressed as

$$ \mathcal{F}_2 \left[ \begin{matrix} a \\ c_1, c_2 \end{matrix} \middle| x, y \right] = \sum_{m,n=0}^{\infty} \Gamma \left[ \begin{matrix} a+m+n, b_1+m, b_2+n \\ c_1+m, c_2+n \end{matrix} \right] \frac{x^m y^n}{m!n!}, \quad (A.8) $$

$$ \mathcal{F}_4 \left[ \begin{matrix} a, b \\ c_1, c_2 \end{matrix} \middle| x, y \right] = \sum_{m,n=0}^{\infty} \Gamma \left[ \begin{matrix} a+m+n, b+m+n \\ c_1+m, c_2+n \end{matrix} \right] \frac{x^m y^n}{m!n!}. \quad (A.9) $$

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

$$ F_4 \left[ \begin{matrix} \frac{a}{2}, \frac{a+1}{2} \\ b_1 + \frac{1}{2}, b_2 + \frac{1}{2} \end{matrix} \middle| x^2, y^2 \right] = (1+x+y)^{-a} F_2 \left[ \begin{matrix} a \\ 2b_1, 2b_2 \end{matrix} \middle| \frac{2x}{x+y+1}, \frac{2y}{x+y+1} \right], \quad (A.10) $$

and the Appell $F_2$ can be expanded as summation of Gauss hypergeometric function $_2F_1$ as

$$ F_2 \left[ \begin{matrix} a \\ c_1, c_2 \end{matrix} \middle| x, y \right] = \sum_{m=0}^{\infty} \frac{(a)_{m} (b_1)_m}{(c_1)_m} \frac{x^m}{m!} {}_2F_1 \left[ \begin{matrix} a+m, b_2 \\ c_2 \end{matrix} \middle| y \right]. \quad (A.11) $$

Finally, the $F_C$-type of Lauricella's hypergeometric functions of $m$ variables $x_1, \dots, x_m$ appears in appendix D, which is defined by [191],

$$ F_C \left[ \begin{matrix} a, b \\ c_1, \dots, c_m \end{matrix} \middle| x_1, \dots, x_m \right] = \sum_{n_1, \dots, n_m = 0}^{\infty} \frac{(a)_{n_1+\dots+n_m} (b)_{n_1+\dots+n_m}}{(c_1)_{n_1} \dots (c_m)_{n_m} n_1! \dots n_m!} x_1^{n_1} \dots x_m^{n_m}, \quad (A.12) $$
with convergent radius

$$ \sqrt{|x_1|} + \cdots + \sqrt{|x_m|} < 1. \qquad (A.13) $$

Note that Lauricella's $F_C$ function is reduced to Appell's and Gauss hypergeometric functions $F_4$ and $_2F_1$ for $m = 2$ and $m = 1$, respectively.

## 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.47), (2.49), and (2.50). 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.41) by changing the variables

$$ z_1 \equiv -k_{124}\tau_1, \quad z_2 \equiv -k_3\tau_2, \qquad (B.1) $$

and introducing

$$ D_{ab}^{\alpha}(k_3; \tau_1, \tau_2) = k_3^{-3} \hat{D}_{ab}^{\alpha} \left( \frac{R}{2-R} z_1, z_2 \right), \qquad (B.2) $$

where

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

Under these change, eq. (2.41) can be written as

$$ \mathcal{I}_{ab,\alpha}^{p_1 p_2}(R) = H^{-2} \left(\frac{R}{2-R}\right)^{1+p_1} (-ab) \int_0^{\infty} dz_1 dz_2 z_1^{p_1} z_2^{p_2} e^{-iaz_1 - ibz_2} \hat{D}_{ab}^{\alpha} \left(\frac{R}{2-R} z_1, z_2 \right). \qquad (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

$$ \mathcal{I}_{\pm\mp,\alpha}^{p_1 p_2}(R) = \frac{e^{\mp i \frac{\pi}{2} \bar{p}_{12}}}{2^{6+p_{12}-2i\mu_{\alpha}} \pi^{\frac{1}{2}}} \Gamma\left[\frac{5}{2}+p_2+i\mu_{\alpha}, \frac{5}{2}+p_2-i\mu_{\alpha}, \frac{5}{2}+p_1-i\mu_{\alpha}, i\mu_{\alpha}\right] \\ \times R^{\frac{5}{2}+p_1-i\mu_{\alpha}} {}_2F_1\left[\frac{1}{2}-i\mu_{\alpha}, \frac{5}{2}+p_1-i\mu_{\alpha} \middle| R\right] + (\mu_{\alpha} \rightarrow -\mu_{\alpha}), \qquad (B.5) $$

$$ \mathcal{I}_{\pm\pm,\alpha}^{p_1 p_2}(R) = \left\{ \frac{\pm i e^{\mp i \frac{\pi}{2} p_{12}} e^{\mp \pi \mu_{\alpha}}}{2^{6+p_{12}-2i\mu_{\alpha}} \pi^{\frac{1}{2}}} \Gamma\left[\frac{5}{2}+p_2+i\mu_{\alpha}, \frac{5}{2}+p_2-i\mu_{\alpha}, \frac{5}{2}+p_1-i\mu_{\alpha}, i\mu_{\alpha}\right] \right. \\ \times R^{\frac{5}{2}+p_1-i\mu_{\alpha}} {}_2F_1\left[\frac{1}{2}-i\mu_{\alpha}, \frac{5}{2}+p_1-i\mu_{\alpha} \middle| R\right] + (\mu_{\alpha} \rightarrow -\mu_{\alpha}) \\ \left. + \frac{e^{\mp i \frac{\pi}{2} p_{12}} \Gamma(p_{12}+5)}{\mu_{\alpha}^2 + \left(\frac{5}{2}+p_2\right)^2} \left(\frac{R}{2}\right)^{5+p_{12}} {}_3F_2\left[\frac{1, 5+p_{12}, 3+p_2}{2}+i\mu_{\alpha}, \frac{7}{2}+p_2+i\mu_{\alpha} \middle| R\right], \right\} \qquad (B.6) $$

where $p_{12} \equiv p_1 + p_2$ and $\bar{p}_{12} \equiv p_1 - p_2$.
### C Boundary conditions of seed integrals

In the main text, in section 3.2, we imposed the soft limit $(u, v \to 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,

$$H_{i\mu}^{(j)}(-k\tau) = \int_{-i\infty}^{i\infty} \frac{ds}{2\pi i} \frac{(-k\tau/2)^{-2s}}{\pi} e^{(-1)^{j+1}(2s-i\mu-1)i\frac{\pi}{2}} \Gamma\left[s-\frac{i\mu}{2}, s+\frac{i\mu}{2}\right], \quad (j=1,2) \quad (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

$$\begin{aligned} D_{\pm\mp}^{\alpha}(k_1; \tau_1, \tau_2) &= \frac{H^2 \pi e^{-\pi \mu_{\alpha}}}{4} (\tau_1 \tau_2)^{3/2} H_{\mp i \mu_{\alpha}}^{(2,1)}(-k_1 \tau_1) H_{\pm i \mu_{\alpha}}^{(1,2)}(-k_1 \tau_2) \\ &= \frac{H^2}{4\pi} \int_{-i\infty}^{i\infty} \frac{ds_1}{2\pi i} \frac{ds_2}{2\pi i} e^{\mp i\pi(s_1-s_2)} \left(\frac{k_1}{2}\right)^{-2s_{12}} (-\tau_1)^{-2s_1+3/2} (-\tau_2)^{-2s_2+3/2} \\ & \qquad \times \Gamma\left[s_1 - \frac{i\mu_{\alpha}}{2}, s_1 + \frac{i\mu_{\alpha}}{2}, s_2 - \frac{i\mu_{\alpha}}{2}, s_2 + \frac{i\mu_{\alpha}}{2}\right], \end{aligned} \quad (C.2)$$

and similarly for $D_{\pm\mp}^{\beta}(k_3; \tau_2, \tau_3)$. For time-ordered propagators, it has been observed that a more effective way involves dividing them into factorised and time-ordered components through

$$D_{\pm\pm}^{\alpha}(k_1; \tau_1, \tau_2) = D_{\pm\mp}^{\alpha}(k_1; \tau_1, \tau_2) + [D_{\pm\pm}^{\alpha}(k_1; \tau_1, \tau_2) - D_{\pm\mp}^{\alpha}(k_1; \tau_1, \tau_2)] \theta(\tau_1 - \tau_2), \quad (C.3)$$

$$D_{\pm\pm}^{\beta}(k_3; \tau_2, \tau_3) = D_{\pm\mp}^{\beta}(k_3; \tau_2, \tau_3) + [D_{\pm\mp}^{\beta}(k_3; \tau_2, \tau_3) - D_{\pm\pm}^{\beta}(k_3; \tau_2, \tau_3)] \theta(\tau_3 - \tau_2). \quad (C.4)$$

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

$$\mathcal{I}_{\pm\mp\pm}, \quad (C.5)$$

$$\mathcal{I}_{\pm\pm\mp} = \mathcal{I}_{\pm\pm\mp,F} + \mathcal{I}_{\pm\pm\mp,T}, \quad (C.6)$$

$$\mathcal{I}_{\pm\mp\mp} = \mathcal{I}_{\pm\mp\mp,F} + \mathcal{I}_{\pm\mp\mp,T}, \quad (C.7)$$

$$\mathcal{I}_{\pm\pm\pm} = \mathcal{I}_{\pm\pm\pm,FF} + \mathcal{I}_{\pm\pm\pm,TF} + \mathcal{I}_{\pm\pm\pm,FT} + \mathcal{I}_{\pm\pm\pm,TT}, \quad (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

$$\begin{aligned} \mathcal{I}_{\pm\mp\pm} = \pm i H^{-4} k_{24}^{9+p123} & \int_{-\infty}^{0} d\tau_1 d\tau_2 d\tau_3 (-\tau_1)^{p_1} (-\tau_2)^{p_2} (-\tau_3)^{p_3} e^{\pm i k_1 \tau_1 \mp i k_{24} \tau_2 \pm i k_3 \tau_3} \\ & \times D_{\pm\mp}^{\alpha}(k_1; \tau_1, \tau_2) D_{\pm\pm}^{\beta}(k_3; \tau_2, \tau_3), \end{aligned} \quad (C.9)$$
$$ \mathcal{I}_{\pm\pm\mp,F} = \pm iH^{-4}k_{24}^{9+p_{123}} \int_{-\infty}^{0} d\tau_1 d\tau_2 d\tau_3 (-\tau_1)^{p_1} (-\tau_2)^{p_2} (-\tau_3)^{p_3} e^{\pm ik_1\tau_1 \pm ik_{24}\tau_2 \mp ik_3\tau_3} \\ \times D_{\pm\mp}^{\alpha}(k_1; \tau_1, \tau_2) D_{\pm\mp}^{\beta}(k_3; \tau_2, \tau_3), \qquad (C.10) $$

$$ \mathcal{I}_{\pm\mp\mp,F} = \mp iH^{-4}k_{24}^{9+p_{123}} \int_{-\infty}^{0} d\tau_1 d\tau_2 d\tau_3 (-\tau_1)^{p_1} (-\tau_2)^{p_2} (-\tau_3)^{p_3} e^{\pm ik_1\tau_1 \mp ik_{24}\tau_2 \mp ik_3\tau_3} \\ \times D_{\pm\mp}^{\alpha}(k_1; \tau_1, \tau_2) D_{\pm\mp}^{\beta}(k_3; \tau_2, \tau_3), \qquad (C.11) $$

$$ \mathcal{I}_{\pm\pm\pm,FF} = \mp iH^{-4}k_{24}^{9+p_{123}} \int_{-\infty}^{0} d\tau_1 d\tau_2 d\tau_3 (-\tau_1)^{p_1} (-\tau_2)^{p_2} (-\tau_3)^{p_3} e^{\pm ik_1\tau_1 \pm ik_{24}\tau_2 \pm ik_3\tau_3} \\ \times D_{\pm\mp}^{\alpha}(k_1; \tau_1, \tau_2) D_{\mp\pm}^{\beta}(k_3; \tau_2, \tau_3). \qquad (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 get 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

$$ \mathcal{I}_{\pm\pm\pm,\alpha\beta}^{p_1 p_2 p_3} = \frac{e^{\mp 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[ \frac{5}{2} + p_1 - i\mu_\alpha, \frac{5}{2} + p_1 + i\mu_\alpha, \frac{5}{2} + p_3 - i\mu_\beta, \frac{5}{2} + p_3 + i\mu_\beta \right] \\ 3 + p_1, 3 + p_3 \\ \times \sum_{\mathbf{a},\mathbf{b}=\pm} \sum_{n_2, n_3=0}^{\infty} \left(\frac{u}{2}\right)^{2n_2 - i\mu_\alpha} \left(\frac{v}{2}\right)^{2n_3 - i\mu_\beta} \frac{(-1)^{n_2}}{n_2! n_3!} \\ \times \Gamma[-n_2 + i\mathbf{a}\mu_\alpha, -n_3 + i\mathbf{b}\mu_\beta, p_2 + 4 + 2n_2 - i(\mathbf{a}\mu_\alpha + \mathbf{b}\mu_\beta)], \qquad (C.13) $$

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,

$$ \Gamma[z, 1-z] = \frac{\pi}{\sin \pi z}, \quad \Gamma \left[ \frac{2z}{z, z + \frac{1}{2}} \right] = \frac{2^{2z-1}}{\pi^{1/2}}, \qquad (C.14) $$

one may summarize

$$ \lim_{u,v \ll 1} \mathcal{I}_{\pm\pm\pm,\alpha\beta}^{p_1 p_2 p_3} = - \sum_{\mathbf{a},\mathbf{b}=\pm} e^{\mp i \frac{\pi}{2} (p_{13} - p_2)} \text{csch}(\pi \mathbf{a} \mu_\alpha) \text{csch}(\pi \mathbf{b} \mu_\beta) \tilde{\Gamma}(p_1, p_2, p_3, \mu_\alpha, \mu_\beta) \\ \times u^{-\frac{5}{2}-p_1-i\mathbf{a}\mu_\alpha} v^{-\frac{5}{2}-p_3-i\mathbf{b}\mu_\beta} \Gamma \left[ \frac{4+p_2-i(\mathbf{a}\mu_\alpha+\mathbf{b}\mu_\beta)}{2}, \frac{5+p_2-i(\mathbf{a}\mu_\alpha+\mathbf{b}\mu_\beta)}{2} \right], \qquad (C.15) $$

with

$$ \tilde{\Gamma}(p_1, p_2, p_3, \mu_\alpha, \mu_\beta) = \frac{\pi^{\frac{1}{2}}}{2^{4+p_{13}-p_2}} \Gamma \left[ \frac{5}{2} + p_1 - i\mu_\alpha, \frac{5}{2} + p_1 + i\mu_\alpha, \frac{5}{2} + p_3 - i\mu_\beta, \frac{5}{2} + p_3 + i\mu_\beta \right] \\ 3 + p_1, 3 + p_3 \qquad (C.16) $$
A similar procedure can be applied to both the partially factorised (partially nested) and fully nested seed integrals. Under double soft limits, they yield:

$$ \lim_{u,v \ll 1} \mathcal{I}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3} = \mp i \sum_{a,b=\pm} e^{\mp i \frac{\pi}{2} (p_{12} - p_3)} \text{csch}(\pi a \mu_{\alpha}) \text{csch}(\pi b \mu_{\beta}) e^{\mp \pi a \mu_{\alpha}} \tilde{\Gamma}(p_1, p_2, p_3, \mu_{\alpha}, \mu_{\beta}) \\ \times u^{-\frac{5}{2}-p_1-ia\mu_{\alpha}} v^{-\frac{5}{2}-p_3-ib\mu_{\beta}} \Gamma\left[ \frac{4+p_2-ia\mu_{\alpha}+b\mu_{\beta}}{2}, \frac{5+p_2-ia\mu_{\alpha}+b\mu_{\beta}}{2} \middle| 1-ia\mu_{\alpha}, 1-ib\mu_{\beta} \right], \quad (C.17) $$

$$ \lim_{u,v \ll 1} \mathcal{I}_{\pm\pm,\alpha\beta}^{p_1 p_2 p_3} = \sum_{a,b=\pm} e^{\mp i \frac{\pi}{2} p_{123}} \text{csch}(\pi a \mu_{\alpha}) \text{csch}(\pi b \mu_{\beta}) e^{\mp \pi (a \mu_{\alpha} + b \mu_{\beta})} \tilde{\Gamma}(p_1, p_2, p_3, \mu_{\alpha}, \mu_{\beta}) \\ \times u^{-\frac{5}{2}-p_1-ia\mu_{\alpha}} v^{-\frac{5}{2}-p_3-ib\mu_{\beta}} \Gamma\left[ \frac{4+p_2-ia\mu_{\alpha}+b\mu_{\beta}}{2}, \frac{5+p_2-ia\mu_{\alpha}+b\mu_{\beta}}{2} \middle| 1-ia\mu_{\alpha}, 1-ib\mu_{\beta} \right]. \quad (C.18) $$

The results for $\mathcal{I}_{\pm\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,\alpha\beta}^{p_1 p_2 p_3}$.

## 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

$$ \mathcal{I}_{abcd, \alpha \beta \gamma}^{p_1 p_2 p_3 p_4} = H^{-6} k_4^{13+p_{1234}} (abcd) \int_{-\infty}^{0} \prod_{i=1}^{4} d\tau_i (-\tau_i)^{p_i} e^{ia k_1 \tau_1 + ib k_2 \tau_2 + ic k_3 \tau_3 + id k_4 \tau_4} \\ \times D_{\alpha d}^{\alpha} (k_1; \tau_1, \tau_4) D_{\beta d}^{\beta} (k_2; \tau_2, \tau_4) D_{\gamma d}^{\gamma} (k_3; \tau_3, \tau_4). \quad (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 \to 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.

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

$$ -k_i \tau_i = z_i, \quad (i = 1, 2, 3, 4), \quad (D.2) $$
<figure><img src="image_16.png" /><figcaption>**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.</figcaption></figure>

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

$$ \begin{aligned} \mathcal{I}_{abcd}^{p_1 p_2 p_3 p_4} &= \frac{H^{-6}(abcd)}{w_1^{4+p_1} w_2^{4+p_2} w_3^{4+p_3}} \int_0^\infty \prod_{i=1}^4 dz_i z_i^{p_i} e^{-iaz_1 - ibz_2 - icz_3 - idz_4} \hat{D}_{ad}(z_1, w_1 z_4) \hat{D}_{bd}(z_2, w_2 z_4) \hat{D}_{cd}(z_3, w_3 z_4), \end{aligned} \quad (D.3) $$

$$ = \frac{1}{w_1^{4+p_1} w_2^{4+p_2} w_3^{4+p_3}} \hat{\mathcal{I}}_{abcd}^{p_1 p_2 p_3 p_4} (w_1, w_2, w_3), \quad (D.4) $$

where

$$ w_i \equiv \frac{k_i}{k_4}, \quad (i = 1, 2, 3). \quad (D.5) $$

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

$$ \mathcal{D}_{w_i}^{(3)} \hat{\mathcal{I}}_{\pm\pm\mp\mp}^{p_1 p_2 p_3 p_4} = 0, \quad (i = 1, 2, 3), \quad (D.6) $$

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

$$ \begin{aligned} \mathcal{D}_{w_i}^{(3)} &= 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) \\ &= w_i^2 (1 - w_i^2) \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 (1 + (p_4 + 2)w_i^2) \partial_{w_i} \\ &\quad - 2(p_4 + 2)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, \end{aligned} \quad (D.7) $$

where $\mu_i$ is the mass parameter of the massive field associated with the momentum $k_i$. Remarkably, this can be analytically solved by

$$ \hat{\mathcal{I}}_{\pm\pm\mp\mp}^{p_1 p_2 p_3 p_4} = \sum_{abc=\pm} c_{\pm\pm\mp\mp,abc} w_1^{\frac{3}{2}-ia\mu_1} w_2^{\frac{3}{2}-ib\mu_2} w_3^{\frac{3}{2}-ic\mu_3} \\ \times F_C \left[ \frac{2p_4+11-2i(a\mu_1+b\mu_2+c\mu_3)}{4}, \frac{2p_4+13-2i(a\mu_1+b\mu_2+c\mu_3)}{4} \middle| w_1^2, w_2^2, w_3^2 \right], \quad (D.8) $$