# The Exact and Approximate Tales of Boost-Breaking Cosmological Correlators

Zhehan Qin, Sébastien Renaux-Petel, Xi Tong, Denis Werth, Yuhang Zhu

► To cite this version:

Zhehan Qin, Sébastien Renaux-Petel, Xi Tong, Denis Werth, Yuhang Zhu. The Exact and Approximate Tales of Boost-Breaking Cosmological Correlators. Journal of Cosmology and Astroparticle Physics, 2025, 09, pp.058. 10.1088/1475-7516/2025/09/058. hal-05116992

HAL Id: hal-05116992

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

Submitted on 31 Oct 2025

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>
PAPER • OPEN ACCESS

# The exact and approximate tales of boost-breaking cosmological correlators

To cite this article: Zhehan Qin et al JCAP09(2025)058

View the [article online](#) for updates and enhancements.

You may also like

- [Building a boostless bootstrap for the bispectrum](#)
  Enrico Pajer

- [Shapes of the cosmological low-speed collider](#)
  Sadra Jazayeri, Sébastien Renaux-Petel and Denis Werth

- [Fingerprints of a non-inflationary universe from massive fields](#)
  Jerome Quintin, Xingang Chen and Reza Ebadi
# The exact and approximate tales of boost-breaking cosmological correlators

Zhehan Qin<sup>DID</sup>,<sup>a,b</sup> Sébastien Renaux-Petel<sup>ID</sup>,<sup>c</sup> Xi Tong<sup>ID</sup>,<sup>b</sup> Denis Werth<sup>IDc</sup>
and Yuhang Zhu<sup>IDd</sup>

**ABSTRACT:** Cosmological correlators offer a remarkable window into the high-energy physics governing Universe's earliest moments, with the tantalising prospect of discovering new particles. However, extracting new physics from these observables requires both precise theoretical predictions of inflationary theories and accurate, analytical templates suitable for data analysis throughout parameter and kinematic spaces. In this paper, we extend the current analytical results by computing the most general boost-breaking seed correlator mediated by the tree-level exchange of a massive spinning particle. We derive the result using two complementary approaches, bootstrapping from boundary differential equations, and direct spectral integration. Both representations are packaged as a single partially resummed series that converges in all physical kinematics. Computing this correlator marks a milestone for carving out the space of all boost-breaking correlators, and therefore completes the tree-level catalogue. We then introduce a general procedure to obtain accurate approximations for cosmological collider signals based on the saddle-point method. This approach allows for a clear physical intuition of various signals hidden in correlators, as the bulk physics is made manifest through the location of these saddles in the complex time plane, which depend on the external kinematics. Evaluating the time integrals at these saddles yields results given as elementary functions that remain valid beyond soft limits and provide intuitive control over both the signal shape and amplitude. We demonstrate the power of this method in both de Sitter-invariant and boost-breaking scenarios, and uncover novel refined waveform and strength dependence for oscillatory signals from massive fields. We provide a complete cosmological collider shape template capturing all boost-breaking effects for upcoming cosmological surveys.

**KEYWORDS:** cosmological perturbation theory, inflation, non-gaussianity, cosmology of theories beyond the SM

ARXIV EPRINT: 2506.01555

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

<footnote><sup>a</sup>Department of Physics, Tsinghua University, Beijing 100084, China</footnote>

<footnote><sup>b</sup>Department of Applied Mathematics and Theoretical Physics,
Cambridge University, Cambridge, CB3 0WA, U.K.</footnote>

<footnote><sup>c</sup>Institut d'astrophysique de Paris, Sorbonne Université, CNRS, 98 bis bd Arago, Paris, France</footnote>

<footnote><sup>d</sup>Cosmology, Gravity and Astroparticle Physics Group,
Center for Theoretical Physics of the Universe, Institute for Basic Science, Daejeon 34126, Korea
E-mail: qzh21@mails.tsinghua.edu.cn, renaux@iap.fr, xt246@cam.ac.uk,
werth@iap.fr, yhzhu@ibs.re.kr</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>1.1</td><td>Exact boost-breaking catalogue</td><td>2</td></tr><tr><td>1.2</td><td>Saddling up for approximate correlators</td><td>3</td></tr><tr><td>1.3</td><td>Summary of main results</td><td>5</td></tr><tr><td>2</td><td>Bootstrapping boost-breaking cosmological correlators</td><td>7</td></tr><tr><td>2.1</td><td>Boost-breaking physics & helical seed correlator</td><td>8</td></tr><tr><td>2.2</td><td>Via bootstrap equation</td><td>10</td></tr><tr><td>2.3</td><td>Via spectral representation</td><td>18</td></tr><tr><td>3</td><td>Classifying non-analyticities</td><td>27</td></tr><tr><td>3.1</td><td>Extracting & approximating non-analyticities</td><td>27</td></tr><tr><td>3.2</td><td>Application to tree-level exchange correlators</td><td>30</td></tr><tr><td>3.3</td><td>Refined cosmological collider templates</td><td>46</td></tr><tr><td>4</td><td>Conclusions</td><td>51</td></tr><tr><td>A</td><td>Mathematical interlude</td><td>52</td></tr><tr><td>A.1</td><td>Riemann surfaces from multi-valued functions</td><td>53</td></tr><tr><td>A.2</td><td>Saddle point approximation of complex integrals</td><td>56</td></tr></table>

# 1 Introduction

The early Universe serves as a unique laboratory for testing the fundamental laws of Nature at energy scales far beyond the reach of terrestrial experiments. In particular, by magnifying quantum fluctuations from microscopic to macroscopic scales, cosmic inflation allows to probe quantum physics through cosmological observations. Among the most promising avenues for uncovering new physics are higher-point cosmological correlators of primordial fluctuations, commonly referred to as primordial non-Gaussianities. These detection channels are a central target of current and next-generation experiments, including observations of the cosmic microwave background and large-scale structure surveys, see [1, 2] for reviews and the references therein.

A close interplay has emerged between theoretical high-energy physics, which aims to compute and classify the space of cosmological correlators, and observational cosmology, which seeks to constrain them with increasing precision. This synergy is embodied in the “cosmological collider program” [3–7], see also [8–40]. Much like in ground-based colliders, where particle interactions are inferred from scattering data, cosmological correlators encode the signatures of particles that were spontaneously produced during inflation and decayed into curvature perturbations. Since the dynamics of the inflationary bulk is ultimately encoded in the spatial patterns of density fluctuations observed at late times, analysing the kinematic dependence of cosmological correlators enables a form of particle spectroscopy —
revealing the masses, spins, dispersion relations and couplings of particles that were active
in the primordial Universe.

To advance this program, it is essential to develop predictions for cosmological correlators
that are *precise*, *efficient*, and *complete*, in the sense of covering wide regions of parameter
space and being valid across all kinematic regimes. In this paper, we address these needs
by completing the catalogue of boost-breaking cosmological correlators and by developing
a systematic framework for constructing approximate solutions. Our approach ultimately
provides refined templates for cosmological collider signals, facilitating their use in data
analysis.

## 1.1 Exact boost-breaking catalogue

Significant progress has been made in our theoretical understanding of cosmological correlators, with exact analytical results now available for a wide range of physical processes. This progress has been driven by the development of powerful computational techniques, including solving boundary differential equations in kinematic space [41–48], the use of Mellin-space representations of massive mode functions exploiting approximate dilatation invariance [24, 44, 49–54], as well as spectral or dispersive methods [55–57]. As a result, we now have access to fully analytical expressions for the shape of some cosmological correlators across all kinematic configurations. Nevertheless, carving out most of the theory space achieving a systematic classification of the full range of possible signatures remains a major challenge.

Building novel computational techniques has been closely guided by symmetries. A key
insight is that, to first order in the mixing with dynamical gravity, inflationary fluctuations can
be well approximated by weakly coupled quantum fields in a fixed (quasi) de Sitter background.
Exploiting the full de Sitter symmetry group — importantly including de Sitter boosts —
are known to completely fix many correlators and enabled the derivation of the first exact,
closed-form solutions for cosmological correlators involving massive particle exchange [6, 40, 57–63].
These developments closely parallel how conformal symmetry constrains correlators in
conformal field theories, and how Poincaré symmetry shapes scattering amplitudes in flat space.
Introducing a mild breaking of boost symmetry allowed these results to be extended to slow-roll
inflationary scenarios, leading to predictions for primordial non-Gaussianities with amplitudes
comparable to the gravitational floor in equilateral configurations, $f_{\text{NL}}^{\text{eq}} = \mathcal{O}(\epsilon, \eta)$, where
$\epsilon, \eta \sim 10^{-2}$ denote the slow-roll parameters [64]. Notably, this predicts a cosmological collider
signal lying in soft limits $k_L/k_S \ll 1$ — where $k_L$ (resp. $k_S$) denotes some long (resp. short)
wavelength momentum — whose amplitude is highly suppressed i.e. $f_{\text{NL}}^{\text{CC}} = \mathcal{O}(e^{-\pi\mu})$, where $\mu$
denotes the mass of the exchanged heavy particle in Hubble units.<fnref n="1" /> However from observations,
we expect primordial fluctuations to be statistically homogeneous, isotropic and approximately
scale-invariant, but not necessarily de Sitter boost-invariant. In fact, we generally expect
the rest frame of the inflaton background to single out a preferred direction of time and
break de Sitter boosts. This distinction is crucial: phenomenologically, the size of primordial
non-Gaussianities is directly related to the amount of symmetry breaking, and rightfully
relaxing certain symmetries not only enlarges the space of inflationary models, but also opens

<footnote n="1">¹This is sometimes called a Boltzmann suppression since the Boltzmann factor corresponding to a de Sitter space of Gibbons-Hawking temperature $T = H/(2\pi)$ is $e^{-m/T} \approx e^{-2\pi\mu}$, where $\mu^2 = m^2/H^2 - 9/4$.</footnote>
up new regions of kinematic space. Since the rolling of the inflaton background automatically breaks boosts, the most natural way to introduce the boost-breaking effect is starting from a de Sitter-covariant theory, and taking (some of) the inflaton fields to the background value, inducing operators that break boosts. While these boost-breaking operators could be treated perturbatively if the theory is weakly coupled, more nontrivial effects arise when they become non-perturbative. In particular, if these operators are quadratic, resumming them will in general modify the dispersion relation of the field content.

To second-order in (spatial) derivatives in the linear theory, all possible boost-breaking effects can be captured by two key parameters: (i) a non-unit speed of sound $c_s \neq 1$, which distinguishes the sound cone from the light cone, and (ii) a non-zero helical chemical potential $\kappa \neq 0$, which arises when coupling matter fields through the rolling inflaton background. These parameters modify the dispersion relation of the massless field and the massive particle, which schematically take the form $\omega_\varphi^2 \sim c_s^2 k_p^2$ and $\omega_\sigma^2 \sim k_p^2 + \kappa k_p + m^2$, with $k_p \equiv k/a$ being the physical (redshifting) momentum, and $m$ the mass of the heavy field.<fnref n="2" /> Taken individually, the effects of these boost-breaking parameters are well understood. On one hand, a reduced sound speed $c_s < 1$ for the curvature perturbation enhances equilateral-type non-Gaussianities $f_{\text{NL}}^{\text{eq}} = \mathcal{O}(c_s^{-2})$ [65–67], and amplifies the cosmological collider signal, $f_{\text{NL}}^{\text{CC}} = \mathcal{O}(e^{-\pi\mu/2})$ in ultra-soft configurations $k_L/k_S \ll c_s$ [11, 26]. Notably, strong boost breaking via subluminal sound speed allows for the novel mildly-soft kinematic region ($k_L/k_S \approx c_s m/H$) to emerge, where a new cosmological low-speed collider resonance appears [26, 43]. Closed-form solutions for the corresponding cosmological correlators were derived in [42, 43]. On the other hand, a non-zero chemical potential is known to enhance particle production [68, 69], see also [14, 17, 18, 24, 30, 44, 70–76]. A direct consequence is the exponential amplification of the cosmological collider signal $f_{\text{NL}}^{\text{CC}} = \mathcal{O}(e^{-\pi(\mu-\kappa)})$ [17]. The chemical potential is sensitive to the helical structure of the exchanged massive spinning field and therefore also leads to parity violation in cosmological correlators [15, 27, 77]. The exact analytical form of helical cosmological correlators have been derived in [44, 45].

Here, we combine the leading boost-breaking effects and compute the four-point seed correlator of external legs having a reduced sound speed, mediated by the tree-level exchange of massive and helical spinning particles. We present two complementary representations of this correlator: one obtained by the bootstrap approach, and the other obtained by extending the recently developed spectral method. Both techniques yield partially resummed single-series solutions that converge across all kinematic configurations. This result completes the tree-level catalogue of boost-breaking cosmological correlators.<fnref n="3" />

## 1.2 Saddling up for approximate correlators

Even at tree level, cosmological correlators involving massive particle exchange are notoriously complex. Their exact expressions typically involve hypergeometric functions and infinite series that must be analytically continued to cover the full kinematic domain. Despite

<footnote n="2">²Without loss of generality for our purposes, we rescaled the sound speed of the heavy field to unity. However, note that a non-unit sound speed necessarily comes with cubic interactions with the massless Goldstone field, with important observational consequences in general, see [28].</footnote>

<footnote n="3">³While we mainly focus on the trispectrum, it is convenient to derive the bispectrum, which is more observationally relevant, by taking suitable soft limits.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Illustration of the simplified, refined, and exact cosmological collider waveform templates (see table 1) as function of the momentum ratio $k_L/k_S$ in the soft limit, varying the mass parameter.</figcaption></figure>

their generality and precision, these analytical solutions are nevertheless complicated to manipulate and the physical picture often remains concealed by obscure mathematics. For instance, it remains conceptually unclear why the amplitudes of cosmological collider signals are always proportional to certain fractional powers of the characteristic Boltzmann factor, i.e. $f_{\text{NL}}^{\text{CC}} = \mathcal{O}(e^{-\pi\mu/n})$, with $n = 1, 2$. More critically, these expressions pose serious challenges for practical applications. Accurate numerical evaluation is computationally expensive, limiting their use in data analysis pipelines. Two main obstacles hinder the detection and constraint of cosmological collider signals: (i) the *exact* waveforms for the oscillatory features from massive particle production involve complicated special functions that are slow to evaluate compared to elementary functions; (ii) the portion of the signal that is observationally accessible lies in the mildly-soft regime — as realistic observations cannot probe arbitrarily soft modes due to cosmic variance — where it is dwarfed by large equilateral contaminations from curvature perturbation self-interactions or more generally effective field theory contributions. These “background signals” must be subtracted with high precision, yet they too are typically expressed as slowly converging infinite series of special functions.

Fortunately, the signature of massive particles in cosmological correlators are among the most transparent. Because particles of different de Broglie wavelengths born and decay at different times during inflation, their dynamics becomes imprinted in specific kinematic configurations of late-time correlators. In particular, the soft limit isolates a regime where the massive mode can be expanded at late times: it oscillates with both positive and negative frequency components due to particle production while being redshifted by the expanding background. This clear physical picture enables the construction of *simplified* waveform templates for cosmological collider signals, as illustrated in figure 1. These signals typically oscillate as a function of the logarithm of the momentum ratio, with a frequency set by the particle's mass: $(k_L/k_S)^{+i\mu} + \text{c.c.} \sim \cos[\mu \log(k_L/k_S)]$, providing a striking observational signature of heavy fields during inflation. Yet, such a simple template is only valid in the ultra-soft limit due to the limitation of late-time approximation of the mode functions. As
<table><caption>**Table 1.** Summary of the key properties — including the template schematic form, kinematic regime of validity, computational cost, and physical intuition — for the three levels of complexity of templates describing the de Sitter invariant cosmological collider signal in the soft limit of cosmological correlators. The kinematic dependence of each template schematic form is expressed in terms of $k_L/k_S$, where $k_L$ represents a long-wavelength mode and $k_S$ a short-wavelength mode. In this work, we propose complete boost-breaking versions of these templates.</caption><thead><tr><th></th><th>Simplified</th><th>Refined</th><th>Exact</th></tr></thead><tbody><tr><th>Template</th><td>cos [μ log(k<sub>L</sub>/k<sub>S</sub>)]</td><td>cos [μ arccosh (k<sub>S</sub>/k<sub>L</sub>)]</td><td>2F<sub>1</sub>[1/2 - iμ/2 + iμ; 1 - k<sub>S</sub>/k<sub>L</sub>/2]</td></tr><tr><th>Method</th><td>late-time expansion</td><td>saddle point</td><td>exact evaluation</td></tr><tr><th>Validity</th><td>k<sub>L</sub>/k<sub>S</sub> &lt;&lt; 1</td><td>k<sub>L</sub>/k<sub>S</sub> &asymp; 1</td><td>k<sub>L</sub>/k<sub>S</sub> &le; 1</td></tr><tr><th>Efficiency</th><td>绿色笑 face</td><td>绿色笑 face</td><td>红色哭 face</td></tr><tr><th>Insight</th><td>绿色笑 face</td><td>绿色笑 face</td><td>红色哭 face</td></tr></tbody></table>

one moves away from the ultra-soft limit, significant dephasing error starts to appear, as shown in figure 1.

In search for a clearer bulk physics interpretation and better boundary templates, in this paper, we shall introduce the saddle point method as a powerful tool to classify and gain physical insight into various cosmological collider signals. We show that the production and decay moments correspond to resonant saddles in the time integrals of cosmological correlators. This means that most of the signal is captured by evaluating the time integrals on these saddle points, yielding an accurate expression for the waveform. This approach relies solely on analysing the location of saddle points in the complex time plane, providing a natural explanation for the behaviour of these signals by tracking the motion of these saddles driven by external kinematics. Surprisingly, we find the characteristic rational powers of the Boltzmann factor for cosmological collider signals are solely determined by the arc angle $\theta_n$ touring around the zero of time to reach the complex saddle point, i.e. $f_{\text{NL}}^{\text{CC}} = \mathcal{O}(e^{-\theta_n \mu})$, $\theta_n = \pi/n$, with $n = 1, 2$. By evaluating the time integrals using the saddle point method, we obtain approximate solutions for a range of correlators, expressed entirely in terms of elementary functions. These simplified expressions are not only easy to implement but also offer practical *refined* waveform templates for observational searches. Comparison with other templates is summarised in table 1 and illustrated in figure 1. The novel templates are valid across almost all kinematic and parameter space, reaching a perfect balance between precision, efficiency and completeness.

## 1.3 Summary of main results

For the reader's convenience, our main results are summarised below:
* Exact tale: we provide the exact closed-form solutions for the helical boost-breaking seed correlator mediated by the tree-level exchange of massive and spinning particles. This takes into account the chemical potential effect for the massive mode, and a reduced sound speed for external legs. The correlator is derived using two complementary approaches — the bootstrap and the spectral methods — and is expressed as a single, partially resummed series that converges rapidly across the entire kinematic domain. Both representations are complementary in terms of computational cost and convergence.

* Approximate tale: we introduce a novel general method for constructing approximate solutions to cosmological collider signals, detailed in section 3.1. Based on the resonant production/decay picture of bulk evolution, this approach combines a Wentzel-Kramers-Brillouin (WKB) approximation for the massive internal mode with a saddle-point evaluation of the time integrals, providing an efficient and physically transparent way to capture and classify the non-analytic features of correlators associated with particle production.

* Phenomenology: using the saddle-point approach, we uncover several key phenomenological features arising from boost symmetry breaking, which we summarize below:

- We find a novel cosmological collider waveform template that extends up to equilateral-type configurations. At the level of the bispectrum dimensionless shape function, it takes the form

$$S^{CC} \sim \cos[\mu \operatorname{arccosh}(k_S/k_L) + \delta]. \quad (1.1)$$

- In the case of a reduced sound speed, we find the following novel interpolated expression for the amplitude of the cosmological collider signal

$$f_{\mathrm{NL}}^{\mathrm{CC}} \sim \exp [-\pi \mu / 2 - \mu \arcsin c_s], \quad (1.2)$$

that interpolates between the usual de Sitter invariant case ($e^{-\pi\mu}$), and the ultra-reduced sound speed ($e^{-\pi\mu/2}$).

- In the presence of a non-zero chemical potential, we highlight a novel signature dubbed transient cosmological collider signal, which takes the form

$$S^{\mathrm{CC}} \sim \cos [c_s \kappa k_L / k_S + \delta], \quad (1.3)$$

oscillating linearly in the momentum ratio with a frequency fixed by the sound speed and the chemical potential.

- We provide a complete refined template for the bispectrum shape capturing all boost-breaking effects in (3.82).

**Outline.** The outline of the paper is as follows: in section 2, we start by briefly reviewing the physics of broken boosts, and introduce the boost-breaking helical seed correlator (section 2.1). We then compute this correlator using two different approaches: (i) the bootstrap method (section 2.2), and (ii) the spectral approach (section 2.3). Both results are packaged into a
single partially resummed series that converges in all physical kinematic configurations. Even-
tually, we benchmark both representations in terms of convergence rate, computational cost,
and analytic continuation. In section 3, we introduce the saddle point method supplemented
by the WKB expansion to obtain approximate yet remarkably accurate expressions for the
non-analytic signals hidden in soft limits of cosmological correlators (section 3.1). We then
apply this method to a number of concrete cases with increasing complexity (section 3.2), and
find refined simple templates for cosmological collider signals. Finally, we state our conclusions
and discuss future directions in section 4. Appendix A provides mathematical details about
multi-valued complex functions and Riemann surfaces, along with a pedagogical introduction
to the saddle point method for complex integrals. It includes detailed computations, whose
results are referenced in the main text, ensuring the paper remains self-contained.

**Notation and conventions.** We use the mostly-plus signature for the metric (–, +, +, +). Throughout most of the paper, we work in the Poincaré patch of de Sitter space in cosmic time *t* and conformal time *τ* with the metric

$$
ds^2 = -dt^2 + a^2(t)d\boldsymbol{x}^2 = a^2(\tau)(-d\tau^2 + d\boldsymbol{x}^2), \qquad (1.4)
$$

where $a(t) = e^{Ht} = a(\tau) = -(H\tau)^{-1}$ is the scale factor and $H$ is the Hubble parameter. To avoid the cluttering of symbols, we will often set the Hubble scale to unity i.e. $H = 1$. Spatial three-dimensional vectors are written in boldface $\mathbf{k}$ (with $\hat{\mathbf{k}} \equiv \mathbf{k}/|\mathbf{k}|$), and spatial indices are denoted with Latin letters $i, j, \dots$. Removal of the boldface denotes taking the magnitude i.e. $k = |\mathbf{k}|$. Momenta in plain font with multiple labels denote their sum, e.g. $k_{12} \equiv k_1 + k_2$ and $k_{34} \equiv k_3 + k_4$. The cross-ratios of momenta are defined by

$$
u \equiv \frac{s}{c_s k_{12}}, \quad v \equiv \frac{s}{c_s k_{34}}, \qquad (1.5)
$$

and

$$
\tilde{u} \equiv \frac{2s}{c_s k_{12} + s}, \quad \tilde{v} \equiv \frac{2s}{c_s k_{34} + s}, \qquad (1.6)
$$

with $\mathbf{s} \equiv \mathbf{k}_1 + \mathbf{k}_2$ the $\mathbf{s}$-channel momentum and $\mathbf{s}$ its magnitude. We use sans serif letters ($\mathbf{a}, \mathbf{b}, \dots$) for Schwinger-Keldysh indices. A prime on correlators denotes the stripping of the momentum conserving $\delta$-function,

$$
\langle \mathcal{O}(\{\mathbf{k}\}) \rangle \equiv \langle \mathcal{O}(\{\mathbf{k}\}) \rangle' (2\pi)^3 \delta^3 (\sum \mathbf{k}) . \qquad (1.7)
$$

Additional definitions will be introduced when needed in the main text.

2 Bootstrapping boost-breaking cosmological correlators

The central object of our interest is the boost-breaking four-point seed correlator that
corresponds to massive-particle exchange at tree level. In this section, we bootstrap this
seed using two different approaches. First, we solve the boundary differential equation in
external kinematics satisfied by the seed correlator and fix unconstrained coefficients using
physics inputs — in our case, by matching the internal soft limit. Second, we leverage the
recently developed spectral representation method to compute this seed. Both methods yield
a single-series representation that converges in all physical kinematic configurations. We
eventually compare both representations in terms of evaluation speed and convergence.
## 2.1 Boost-breaking physics & helical seed correlator

We begin by introducing the physics of broken boosts. In this work, we allow external fields to have a reduced speed of sound $c_s < 1$, and consider the exchange of a massive (spinning) particle that contains linear terms in its dispersion relation $\omega^2 \supset k_p$, where $k_p \equiv k/a$ is the physical momentum. This is achieved by coupling the massive field to a helical chemical potential, whose physical origin we review here. We end this subsection by introducing the helical seed correlator, which is the most general boost-breaking correlator that contains all possible boost-breaking effects up to the second order in (spatial) derivatives.

**Spinning helical fields.** Massive spinning particles during inflation can be described and classified based on the unbroken (three-dimensional) rotational symmetry group [78], see also [77] for further details. This approach is particularly well-suited for studying the effects of boost-breaking, providing a more general framework than classifications based on the full de Sitter isometry group, as considered in [11]. Massive spinning fields with integer spin $S$ are defined as totally-symmetric and traceless 3-tensors, $\sigma_{i_1...i_S}$ with $i_1, \dots, i_S$ running over all spatial dimensions. This object contains $(S_2^+ - S_2^-) = 2S + 1$ degrees of freedom. The most general quadratic action up to two derivatives is

$$ S = \frac{1}{2S!} \int dt d^3 x a^3 \left[ \begin{aligned} & \dot{\sigma}_{i_1 \dots i_S}^2 - c^2 \frac{(\partial_j \sigma_{i_1 \dots i_S})^2}{a^2} - \delta c^2 \frac{(\partial_j \sigma_{j_1 \dots i_S})^2}{a^2} \\ & - m^2 \sigma_{i_1 \dots i_S}^2 - 2S \xi \epsilon_{ijk} \sigma_{il_2 \dots l_S} \frac{\partial_j}{a} \sigma_{kl_2 \dots l_S} \end{aligned} \right], \quad (2.1) $$

where we have introduced the four free parameters $c, \delta c, m^2$ and $\xi$, which respectively correspond to two sound speeds, the (squared) mass, and the chemical potential. Note that the spatial indices of $\sigma_{i_1...i_S}$ can be raised with $\delta_{ij}$ as this object lives on the flat constant-time hypersurface. The field $\sigma_{i_1...i_S}(t, \mathbf{x})$ can be expanded into its different helicity components in Fourier space

$$ \sigma_{i_1 \dots i_S}(t, \mathbf{x}) = \int \frac{d^3 k}{(2\pi)^3} \sigma_{i_1 \dots i_S}(t, \mathbf{k}) e^{i\mathbf{k} \cdot \mathbf{x}}, \quad \text{with} \quad \sigma_{i_1 \dots i_S}(t, \mathbf{k}) = \sum_{\lambda=-S}^{S} \sigma_{\mathbf{k}}^{\lambda}(t) \varepsilon_{i_1 \dots i_S}^{\lambda}(\hat{\mathbf{k}}), \quad (2.2) $$

where $\sigma_{\mathbf{k}}^{\lambda}(t)$ are the mode functions and $\varepsilon_{i_1 \dots i_S}^{\lambda}(\hat{\mathbf{k}})$ are the polarisation tensors satisfying the following reality of the field in real space condition and the normalisation choice

$$ [\varepsilon_{i_1 \dots i_S}^{\lambda}(\hat{\mathbf{k}})]^* = \varepsilon_{i_1 \dots i_S}^{\lambda}(-\hat{\mathbf{k}}), \quad \sum_{i_1 \dots i_S} \varepsilon_{i_1 \dots i_S}^{\lambda}(\hat{\mathbf{k}}) [\varepsilon_{i_1 \dots i_S}^{\lambda'}(\hat{\mathbf{k}})]^* = S! \delta_{\lambda \lambda'}. \quad (2.3) $$

The longitudinal polarisation tensor $\varepsilon^0(\hat{\mathbf{k}})$ is a function of $\hat{\mathbf{k}}$ only, whereas the (partially) transverse tensors $\varepsilon^\lambda(\hat{\mathbf{k}}, \hat{\mathbf{e}}^\pm)$ with $\lambda \neq 0$ in addition depend on the transverse polarisation vectors

$$ \hat{\mathbf{e}}^{\pm} = \frac{\hat{n} - (\hat{n} \cdot \hat{k})\hat{k} \pm i\hat{k} \times \hat{n}}{\sqrt{2[1 - (\hat{n} \cdot \hat{k})^2]}}, \quad (2.4) $$

with $\hat{n}$ an arbitrary unit reference vector satisfying $\hat{n} \times \hat{k} \neq 0$. The transverse polarisation vectors are constructed such that $\hat{k} \cdot \hat{\mathbf{e}}^\pm = 0$. Note that the fully transverse tensors depend
only on the transverse vectors i.e. $\epsilon^{\pm S}(\hat{k}, \hat{e}^{\pm}) = \epsilon^{\pm S}(\hat{e}^{\pm})$. We refer the readers to [77] for the detailed form of these polarisation tensors.

In the helical basis, the mode functions decouple from each other and the quadratic action in Fourier space becomes

$$S^{(\lambda)} = \frac{1}{2} \sum_{\lambda=-S}^{+S} \int dt a^3 \int \frac{d^3 k}{(2\pi)^3} \left[ (\dot{\sigma}_k^{\lambda})^2 - c_{S,\lambda}^2 \frac{k^2}{a^2} (\sigma_k^{\lambda})^2 - m^2 (\sigma_k^{\lambda})^2 - 2\lambda\xi \frac{k}{a} (\sigma_k^{\lambda})^2 \right], \quad (2.5)$$

where $c_{S,\lambda}^2 = c^2 + \frac{S^2-\lambda^2}{S(2S-1)}\delta c^2$. The boost-breaking physics is now clear. Each helicity mode can propagate with a different sound speed, and the parity-violating chemical potential term affects all modes with $\lambda \neq 0$ in an asymmetric way. Varying this action yields the equation of motion

$$\left( \frac{\partial^2}{\partial \tau^2} - \frac{2}{\tau} \frac{\partial}{\partial \tau} + c_{S,\lambda}^2 k^2 - 2\lambda \kappa \frac{k}{\tau} + \frac{m^2}{H^2 \tau^2} \right) \sigma_k^{\lambda} = 0, \quad (2.6)$$

where we have switched to conformal time $\tau$, defined by $d\tau \equiv dt/a$, and $\kappa \equiv \xi/H$ is the dimensionless chemical potential in Hubble units. Up to an irrelevant overall phase, and imposing Bunch-Davies initial conditions, the solution is

$$\sigma_k^{\lambda}(\tau) = - \frac{H\tau}{\sqrt{2c_{S,\lambda}k}} e^{-\pi\lambda\kappa/2} W_{i\lambda\kappa,i\mu}(2ic_{S,\lambda}k\tau). \quad (2.7)$$

For later convenience, we also define the tilded symbol $\tilde{\sigma}$ by absorbing one scale factor, i.e. $\tilde{\sigma}_k \equiv a^{-1}(\tau)\sigma_k$, so that

$$\tilde{\sigma}_k^{\lambda}(\tau) = \frac{e^{-\pi \lambda \kappa / 2}}{\sqrt{2c_{S,\lambda}k}} W_{i\lambda\kappa,i\mu}(2ic_{S,\lambda}k\tau), \quad (2.8)$$

where $W_{i\lambda\kappa,i\mu}$ is the Whittaker-W function, and $\mu \equiv \sqrt{m^2/H^2 - 9/4}$ is the dimensionless mass of the $\sigma$ field. For the purpose of studying generic boost-breaking effects, we will focus on the transverse modes ($\lambda \neq 0$) in the following. For a fixed spin $S$ and given that transverse helicities $\pm\lambda$ that propagate with the same effective sound speed — as $c_{S,\lambda}^2 = c_{S,-\lambda}^2$ —, we can rescale spatial coordinates to absorb $c_{S,\lambda}$ in the sound speed of the external field $c_s$, i.e. $c_s \to c_s/c_{S,\lambda} := c_s$. Notice that this choice is degenerate with setting $\delta c = 0$ and $c = 1$, as far as the helicity-$\pm\lambda$ modes are concerned. This way, $c_s < 1$ means that the external field propagates slower than the field $\sigma_{i_1...i_S}$, and the opposite for $c_s > 1$, which by no means implies superluminality. We shall see later that the chemical potential exponentially enhances the cosmological collider signals but could lead to dangerous particle overproduction if it is larger than the mass of the particle [68]. Therefore, in what follows, we restrict ourselves to $\kappa \lesssim \mu$ to avoid such an instability.

**Helical seed integral.** The observables of interest are expectation values of two massless fields, the scalar field inflaton $\varphi$ (or equivalently, the curvature perturbation $\zeta$) and the tensor field graviton $\gamma_{ij}$. In many applications involving helical chemical potentials, it is sufficient to compute a single helical seed integral. The corresponding cosmological correlators can then be easily obtained, as we will demonstrate with an explicit example below. Using the Schwinger-Keldysh (SK) diagrammatic rules [79], one can define the helical seed integral as [45]

$$F_{ab}^{(\lambda)} \equiv -ab s \int_{-\infty}^{0} \frac{d\tau'}{\tau'} e^{ia c_s k_{12} \tau'} \int_{-\infty}^{0} \frac{d\tau''}{\tau''} e^{ib c_s k_{34} \tau''} G_{ab}^{\lambda}(s; \tau', \tau''), \quad (2.9)$$
where $s = |\mathbf{k}_1 + \mathbf{k}_2|$ is the magnitude of the exchanged momentum in the s-channel, and a, b = ± denote SK indices. We have added an overall factor s to ensure the seed integral is dimensionless and scale-invariant. The bulk-bulk SK propagators for the helicity-λ component of the massive spinning field, $G_{ab}^λ$, are constructed out of the mode functions (2.7) and the Heaviside Θ-function,

$$ 
\begin{aligned}
G_{-+}^{\lambda}(k; \tau', \tau'') &= \tilde{\sigma}_k^{\lambda}(\tau') \tilde{\sigma}_k^{\lambda*}(\tau''), & G_{+-}^{\lambda}(k; \tau', \tau'') &= G_{-+}^{\lambda*}(k; \tau', \tau''), \\
G_{\pm\pm}^{\lambda}(k; \tau', \tau'') &= G_{\mp\pm}^{\lambda}(k; \tau', \tau'') \Theta(\tau' - \tau'') + G_{\pm\mp}^{\lambda}(k; \tau', \tau'') \Theta(\tau'' - \tau'),
\end{aligned}
\quad (2.10) $$

**Spin-1 example.** Starting from the boost-breaking helical seed integral (2.9), we can derive a wide range of correlators involving the exchange of massive spinning fields with non-unit sound speeds and chemical potentials. Consider the spin-1 case [27] as a particularly interesting example, the four-point correlator of inflaton fluctuations $\varphi$ is generated by the interaction

$$ S_{\text{int}} = \int d\tau d^3x \left( \frac{a^{-1}}{\Lambda^3} \partial_j \varphi' \partial_i \varphi \sigma_i \right). \quad (2.11) $$

The corresponding four-point exchange correlator can be expressed as

$$ 
\begin{aligned}
\langle \varphi_{k_1} \varphi_{k_2} \varphi_{k_3} \varphi_{k_4} \rangle' = & - \frac{1}{16\Lambda^6} \sum_{\lambda=\pm 1} \sum_{a,b=\pm} (\hat{k}_1 \cdot \hat{k}_2)(\hat{k}_3 \cdot \hat{k}_4) [\hat{k}_1 \cdot \varepsilon^{\lambda}(\hat{s})]^* [\hat{k}_3 \cdot \varepsilon^{\lambda}(\hat{s})] \\
& \times \frac{1}{k_1 k_3 s^7} (1 - k_1 \partial_{k_1})(1 - k_3 \partial_{k_3}) [D_{k_{12}}^a]^3 [D_{k_{34}}^b]^3 F_{ab}^{(\lambda)} + 3 \text{perms.} \\
& + (t\text{-channel}) + (u\text{-channel}),
\end{aligned}
\quad (2.12) $$

where $D_k^a \equiv is/(ac_s)\partial_k$ and a prime on a correlator denotes that we have dropped the momentum conserving delta function This interaction comes from the operator $\nabla_\mu \delta g^{00} \delta K_\nu^\mu \sigma^\nu$ within the EFT of inflation framework [67]. The polarisation factor can be further explicitly written as

$$ [\hat{k}_1 \cdot \varepsilon^{\lambda}(\hat{s})]^*[\hat{k}_3 \cdot \varepsilon^{\lambda}(\hat{s})] = \frac{1}{2} [\hat{k}_1 \cdot \hat{k}_3 - (\hat{k}_1 \cdot \hat{s})(\hat{k}_3 \cdot \hat{s}) + i\lambda \hat{s} \cdot (\hat{k}_1 \times \hat{k}_3)] . \quad (2.13) $$

It is clear from the expression of the kinematic factor (2.13) that the parity-even sector of the four-point correlator is purely real whereas the parity-odd sector is purely imaginary. In the case without a chemical potential, the seed function $F_{ab}^{(\lambda)}$ that encodes dynamical information does not depend on helicity and the parity-odd part in the kinematic factor cancels between positive and negative helicities, leaving the result parity-even. However, a non-zero chemical potential creates an imbalance between the positive and negative helicity modes in their dynamics i.e. $F_{ab}^{(\lambda)} \neq F_{ab}^{(-\lambda)}$, thereby offsetting the parity-odd part of the correlator from zero.

## 2.2 Via bootstrap equation

We now compute the helical seed correlator using the bootstrap approach [41–46]. As we will see, this method yields a double-series representation that only converges in a restricted region of the kinematic space, due to a reduced sound speed $c_s < 1$. However, by resumming one layer of the series, we are able to analytically continue the solution beyond its region of validity to the full physical kinematic region.
**Boundary equation & kinematics.** The helical seed correlator satisfies the following boundary differential equations [45]

$$ \left[ \Delta_{\pm, \tilde{u}}^{(\lambda)} + \left( \mu^2 + \frac{1}{4} \right) \right] F_{\pm \mp}^{(\lambda)} = 0, \quad \left[ \Delta_{\pm, \tilde{u}}^{(\lambda)} + \left( \mu^2 + \frac{1}{4} \right) \right] F_{\pm \pm}^{(\lambda)} = \frac{1}{2} \frac{\tilde{u} \tilde{v}}{\tilde{u} + \tilde{v} - \tilde{u} \tilde{v}}, \quad (2.14) $$

with the differential operator for different helicities defined by

$$ \Delta_{\pm, \tilde{u}}^{(\lambda)} \equiv \tilde{u}^2(1 - \tilde{u})\partial_{\tilde{u}}^2 - [(1 \pm i\lambda\kappa)\tilde{u}^2]\partial_{\tilde{u}}, \quad (2.15) $$

and similar equations for $\tilde{u} \leftrightarrow \tilde{v}$. We have defined the following dimensionless kinematic variables

$$ \tilde{u} \equiv \frac{2s}{c_s k_{12} + s}, \quad \tilde{v} \equiv \frac{2s}{c_s k_{34} + s}. \quad (2.16) $$

Note that these variables differ from the ones introduced in [41–43]

$$ u \equiv \frac{s}{c_s k_{12}}, \quad v \equiv \frac{s}{c_s k_{34}}, \quad (2.17) $$

and are simply related to the tilted variables by the map $\tilde{u} = 2u/(1+u)$ and $\tilde{v} = 2v/(1+v)$. The inverse map is given by $u = \tilde{u}/(2-\tilde{u})$ and $v = \tilde{v}/(2-\tilde{v})$. In the case of a unit sound speed, these two variables are always constrained to be less than unity due to triangle inequality, $0 < \tilde{u}, \tilde{v} < 1$. The solution to (2.14) in this unit-sound-speed case has been found in [44]. However, for a reduced sound speed, $c_s < 1$, $\tilde{u}, \tilde{v}$ may exceed unity i.e. $0 < \tilde{u}, \tilde{v} < 2/(c_s+1)$, which renders the original series solution of [44] superficially divergent.

**Singularity structure.** One advantage of expressing the seed integral (2.9) as a solution to the bootstrap equations (2.14) is that we can readily identify all possible singularities of the correlator in the complex plane, even before deriving its explicit expressions. Since the seed function is symmetric under the exchange $\tilde{u} \leftrightarrow \tilde{v}$, it suffices to fix $\tilde{v}$ in the physical region with $0 < \tilde{v} < 2/(c_s+1)$, and consider the $\tilde{u}$-complex plane. Kinematical singularities can only arise in the following two cases:

* The coefficient of the leading differential operator in $\Delta_{\pm, \tilde{u}}^{(\lambda)}$ vanishes or diverges:

$$ \tilde{u}^2 - \tilde{u}^3 = 0 \quad \text{or} \quad \infty \quad \Rightarrow \quad s = 0 \quad \text{or} \quad c_s k_{12} \pm s = 0. \quad (2.18) $$

The internal soft singularity at $s=0$ corresponds to the branching point of the (non-local) cosmological collider signals. The other two singularities at $c_s k_{12} + s$ and $c_s k_{12} - s$ correspond to the partial-energy and folded poles, respectively.<fnref n="4" />

* The source term diverges:

$$ \tilde{u} + \tilde{v} - \tilde{u} \tilde{v} = 0 \quad \Rightarrow \quad k_{12} + k_{34} = 0, \quad (2.19) $$

which maps to the total-energy singularity.

<footnote n="4">⁴Under the symmetrisation $k_{12} \leftrightarrow k_{34}$, we immediately deduce the presence of the other partial energy singularity at $c_s k_{34} + s = 0$ and folded singularity at $c_s k_{34} - s = 0$.</footnote>
Notice that (2.18) and (2.19) are only necessary conditions for a singular point of the seed integral. As we will see below, the folded poles $c_s k_{12} - s = 0$ (and also $c_s k_{34} - s = 0$) will be eliminated by boundary conditions. Surprisingly, although the sound speed modifies the location of the partial energy singularities and folded singularities in a trivial manner ($k_i \to c_s k_i$ for external fields), neither the chemical potential nor the sound speed alters the singularity structure of the correlator at the tree level. This is rather non-trivial, as both the chemical potential and the sound speed break the spacetime symmetry of the theory, and one might be led to expect the disruption of the analytic structure of correlators. It would be interesting to investigate the underlying reasons and to examine whether the same conclusion holds also at the loop level.

**Homogeneous solution.** The homogeneous equation takes the same form as in the unit sound speed case. Here, we will adhere to the derivation in previous works and present only the key results relevant to our discussion. We refer the interested readers to [45] for more technical details. There are two differential operators that are related to each other through the helicity reflection condition:

$$ \Delta_{-, \tilde{u}}^{(\lambda)} = \Delta_{+, \tilde{u}}^{(-\lambda)}. \qquad (2.20) $$

Since the equation is second order, each differential operator possesses two linearly independent solutions such that

$$ \left[ \Delta_{+, \tilde{u}}^{(\lambda)} + \left( \mu^2 + \frac{1}{4} \right) \right] \mathcal{H}_{\pm}(\tilde{u}) = 0. \qquad (2.21) $$

The homogeneous solutions are explicitly given by [44, 45]

$$ \mathcal{H}_{\pm}(\tilde{u}) = \tilde{u}^{\frac{1}{2} \pm i\mu} \Gamma(\frac{1}{2} \pm i\mu) \Gamma(\frac{1}{2} \pm i\lambda\kappa \pm i\mu) {}_2\tilde{F}_1 \left[ \begin{array}{c} \frac{1}{2} \pm i\mu \quad \frac{1}{2} \pm i\lambda\kappa \pm i\mu \\ 1 \pm 2i\mu \end{array} ; \tilde{u} \right], \qquad (2.22) $$

where ${}_2\tilde{F}_1[a \ b \ c; z] \equiv \frac{1}{\Gamma(c)} {}_2F_1[a \ b \ c; z]$ is the regularised hypergeometric function. The hypergeometric functions exhibit a branch cut along $\tilde{u} \ge 1$. In the case where $c_s = 1$, $\tilde{u}$ remains less than unity, ensuring that the branch point (i.e. the folded configuration) is never crossed in the physical region. However, when $c_s < 1$, the variable $\tilde{u}$ can exceed unity in certain kinematic regions, thereby crossing the branch point. Since the evaluation of the hypergeometric function ${}_2F_1$ for $\tilde{u} > 1$ is well-established, the homogeneous solution remains valid even in this regime. In other words, the homogeneous solution can be smoothly continued from the unit sound speed case to the small sound speed case, owing to the known analytic continuation properties of hypergeometric series. For the factorised contribution, the general solution is then expressed as

$$ F_{+-}^{(\lambda)}(\tilde{u}, \tilde{v}) = A_{++}^{+-} \mathcal{H}_{+}(\tilde{u}) \mathcal{H}_{-}^{*}(\tilde{v}) + A_{+-}^{+-} \mathcal{H}_{+}(\tilde{u}) \mathcal{H}_{+}^{*}(\tilde{v}) \\ + A_{+-}^{+-} \mathcal{H}_{-}(\tilde{u}) \mathcal{H}_{-}^{*}(\tilde{v}) + A_{--}^{+-} \mathcal{H}_{-}(\tilde{u}) \mathcal{H}_{+}^{*}(\tilde{v}), \qquad (2.23) $$

and $\mathcal{A}$ represents some coefficients that will be determined through the boundary conditions later. We naturally have $F_{-+}^{(\lambda)} = F_{+-}^{(\lambda)*}$. Also for the nested contribution, we can express it as

$$ F_{++}^{(\lambda)}(\tilde{u}, \tilde{v}) = A_{++}^{++} \mathcal{H}_{+}(\tilde{u}) \mathcal{H}_{+}(\tilde{v}) + A_{+-}^{++} \mathcal{H}_{+}(\tilde{u}) \mathcal{H}_{-}(\tilde{v}) \\ + A_{+-}^{++} \mathcal{H}_{-}(\tilde{u}) \mathcal{H}_{+}(\tilde{v}) + A_{--}^{++} \mathcal{H}_{-}(\tilde{u}) \mathcal{H}_{-}(\tilde{v}) \\ + S_{+}^{(\lambda)}(\tilde{u}, \tilde{v}), \qquad (2.24) $$

where $S_{+}^{(\lambda)}$ is a particular solution for the inhomogeneous differential equation which we solve now. As for the other nested seed function, it can be obtained using the relation $F_{--}^{(\lambda)} = F_{++}^{(\lambda)*}$.
**Particular solution.** To find a particular solution, we assume $\tilde{u} < \tilde{v}$ without loss of generality. The result for the other region is obtained by swapping the relevant variables. Note that this procedure is highly non-trivial for higher-point correlators as one needs to analytically continue correlators to other kinematic regions, picking up homogeneous solutions along the way. In this discussion, we focus on $S_+$ as an illustrative example, while the corresponding result for $S_-$ (the particular solution of $F^{(\lambda)}$) can be derived straightforwardly using the conjugate relation $S_- = [S_+]^*$. The source term in eq. (2.14) can be systematically expanded as follows

$$ \frac{\tilde{u}\tilde{v}}{\tilde{u} + \tilde{v} - \tilde{u}\tilde{v}} = \sum_{n=0}^{\infty} \tilde{u}^{n+1} \left(1 - \frac{1}{\tilde{v}}\right)^n = \tilde{u} \sum_{n=0}^{\infty} (-1)^n \left(\frac{\tilde{u}}{\tilde{v}}\right)^n (1 - \tilde{v})^n. \quad (2.25) $$

We propose the following ansatz for the particular solution, assuming a functional form that resembles the structure of the source:

$$ S_{+}^{(\lambda)} = \sum_{m,n=0}^{\infty} A_{m,n}^{(\lambda)} (-1)^n \tilde{u}^{m+1} \left(\frac{\tilde{u}}{\tilde{v}}\right)^n (1 - \tilde{v})^n. \quad (2.26) $$

Substituting the proposed ansatz into the differential equation, we obtain a recurrence relation that governs the coefficients

$$ A_{m+1,n}^{(\lambda)} = \frac{(m+n+1)(m+n+1+i\lambda\kappa)}{(m+n+1+\frac{1}{2})^2 + \mu^2} A_{m,n}^{(\lambda)}, \quad A_{0,n}^{(\lambda)} = \frac{1}{2} \frac{1}{(n+\frac{1}{2})^2 + \mu^2}. \quad (2.27) $$

Using this recurrence relation, the coefficients can be iteratively solved. Explicitly, the coefficients are given by [44, 45]

$$ A_{m,n}^{(\lambda)} = \frac{(n+1)_m (n+1+i\lambda\kappa)_m}{2(n+\frac{1}{2}+i\mu)_{m+1}(n+\frac{1}{2}-i\mu)_{m+1}}, \quad (2.28) $$

where $(a)_m \equiv \Gamma(a+n)/\Gamma(a)$ is the Pochhammer symbol. This series expansion is well-behaved only within the region $\tilde{u} < 1$, indicating that the solution derived here appears to be valid exclusively for the case of unit sound speed. In the kinematic region where $\tilde{u} > 1$, the series fails to converge, suggesting that the ansatz employed here is unsuitable for this range. To address this problem, one could try to attack the equations (2.14) with an alternative ansatz that manifestly converges in $\tilde{u} > 1$ (for example, with negative powers of $\tilde{u}$). For the case without the chemical potential, such ansatz has been adopted and solved in [43]. However, in this work, we provide a more straightforward solution, via a direct resummation of (2.26) that smoothly facilitates its analytic continuation to the $\tilde{u} > 1$ regions. To illustrate this process, we observe that the m-layer of the series can be resummed, and the particular solution becomes

$$ S_{\pm}^{(\lambda)} = \frac{\tilde{u}}{2} \sum_{n=0}^{\infty} \frac{(-1)^n}{(n+\frac{1}{2})^2 + \mu^2} \left(\frac{\tilde{u}}{\tilde{v}}\right)^n (1-\tilde{v})^n {}_3F_2\left[1, n+1, n+1 \pm i\lambda\kappa; \frac{3}{2}+n-i\mu, \frac{3}{2}+n+i\mu; \tilde{u}\right]. \quad (2.29) $$

Here, we have summed the $\tilde{u}$-series into a generalised hypergeometric function. This allows us to evaluate the series even for $\tilde{u} > 1$. Given that we have assumed $\tilde{u} < \tilde{v}$ and by definition $0 < \tilde{v} < 2$, it can be explicitly verified that this resummed single-layer series converges throughout the entire kinematic region under consideration. One aspect we have not yet
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The (rescaled) boost-breaking helical seed correlator $F^{(\lambda)}(u, v)$, as function of the dimensionless kinematic variable $u$ with fixed $v = 0.5$ and at fixed mass $\mu = 3$ for $\kappa = 0, \kappa = 0.5$ and $\lambda = -1$ (left panel) and $\lambda = +1$ (right panel). We include terms in the series from $n = 1, \dots, 10$. We show both the full correlator (solid line), and the background series contribution (dashed line) which is almost identical for all parameters. We considered kinematic variables inside the unit circle corresponding to a unit sound speed $c_s = 1$.</figcaption></figure>

addressed is the branch cut of the hypergeometric series, where $\tilde{u} > 1$ along the real axis. Taking $F_{++}^{(\lambda)}$ as an example, the appropriate branch involves selecting the lower-half plane prescription, specifically $\tilde{u} = \tilde{u} - i\epsilon$ with the infinitesimal $\epsilon > 0$ (vice versa, for $F_{--}^{(\lambda)}$ the appropriate choice, consistent with the complex conjugate relation, is to take $\tilde{u} = \tilde{u} + i\epsilon$). This branch choice can be justified by analysing the corresponding bulk integral, which will be discussed in detail later. As an interesting aside, using the variables $(\tilde{u}, \tilde{v})$ instead of $(u, v)$ allows us to obtain a series representation that converges over a broader kinematic region than originally designed for. Indeed, writing $(\tilde{u}/\tilde{v})^n(1-\tilde{v})^n = \tilde{u}^n(\frac{1-\tilde{v}}{\tilde{v}})^n$, we obtain that the series is convergent for $|\tilde{u}(1-\tilde{v})/\tilde{v}| < 1$. This region is larger than $\tilde{u} < \tilde{v}$ and also extends beyond the unit circle for a reduced sound speed. A natural question is whether one can find another set of well-chosen kinematic variables so that the series converges in the entire physical kinematic domain, which would suggest that (both layers of) the series can be fully resummed into known hypergeometric series. We leave these thoughts for a future work.

**Boundary conditions.** We have nearly completed the calculation of the analytical expressions, except for the unfixed coefficients $A$. Several approaches can be used to impose the boundary conditions, such as requiring the absence of spurious poles or considering the factorised limit as discussed previously. For simplicity, here the coefficients are directly determined by matching the expression in the (hierarchical) soft limit $\tilde{u} \ll \tilde{v} \ll 1$, where the bulk integral can be evaluated easily through the IR expansion of the integrand or using the
Mellin-Barnes transformation. The final expressions are found to be [44, 45]:

$$
\begin{align}
A_{++}^{+-} &= -A_{+-}^{+-} = -A_{-+}^{+-} = A_{--}^{+-} = \frac{-e^{-\pi\lambda\kappa}}{4\sinh^2(2\pi\mu)} [\cosh(2\pi\kappa) + \cosh(2\pi\mu)], \nonumber \\
A_{++}^{++} &= -A_{+-}^{++} = \frac{-i\pi e^{-\pi(\lambda\kappa-\mu)} \cosh[\pi(\mu + \lambda\kappa)]}{2\sinh^2(2\pi\mu)\Gamma(\frac{1}{2}+i\lambda\kappa+i\mu)\Gamma(\frac{1}{2}+i\lambda\kappa-i\mu)}, \tag{2.30} \\
A_{--}^{++} &= -A_{-+}^{++} = \frac{-i\pi e^{-\pi(\lambda\kappa+\mu)} \cosh[\pi(-\mu + \lambda\kappa)]}{2\sinh^2(2\pi\mu)\Gamma(\frac{1}{2}+i\lambda\kappa+i\mu)\Gamma(\frac{1}{2}+i\lambda\kappa-i\mu)}. \nonumber
\end{align}
$$

This concludes the derivation of the helical seed correlator involving two boost-breaking effects: (i) the sound speed, and (ii) the chemical potential. The full exact solution is stored in the Mathematica notebook available on the Github repository linked in ref. [93]. This result is valid for all physical kinematic regions and exhibits rapid convergence. In figure 2, we show the boost-breaking helical seed correlator at a fixed mass $\mu$ but varying the chemical potential $\kappa$ and for both helicities $\lambda = \pm 1$. We note that the chemical potential essentially modifies the signal contribution, boosting the amplitude and slightly changing the phase of the cosmological collider signal, but barely affects the background signal. Also, note that the signal coming from one transverse helicity ($\lambda = -1$) is enhanced compared to the other one. We show in figure 3 the seed correlator for kinematic configurations extending outside the unit circle $1 \le u, v$, which confirms that the found solution is well convergent when the external field has a reduced sound speed. However we have noticed that the convergence rate is slower outside the unit circle than inside. We will come back to this point later when comparing with the spectral representation solution. The solution presents a slight kink around $u = 1$ (from the discontinuity of the first derivative) that disappears swiftly as we dial up the number of terms included $N$, and that is less visible as we also increase the mass $\mu$. The characteristic peak for $\mu = 1$ is the low-speed collider signal [26, 27, 43]. Also notice that the enhanced non-analytic contribution for the negative helicity ($\lambda = -1$) is well visible for $\mu = 3$ and $\kappa = 1$ and with an amplitude that increases as $\kappa$ approaches $\mu$. This is because the amplitude of this cosmological collider signal is given by the Boltzmann factor $\sim e^{-\pi(\mu+\lambda\kappa)}$. We will give more details and thoroughly study this signal in section 3.

**Discussion.** Let us make a few remarks concerning the full solution for the boost-breaking helical massive exchange correlator.

* Bulk cutting rules revisited: bulk cutting rules have been introduced in [80] for tree-level diagrams and subsequently generalised to account for loop-level processes in [51, 52, 81] for extracting non-analytic processes. To briefly revisit these cutting rules, as an illustrative example, consider the integral $F_{++}^{(\lambda)}$. First, assuming that $k_{12} > k_{34}$, by manipulating one of the Heaviside theta functions, $\Theta(\tau' - \tau'') = 1 - \Theta(\tau'' - \tau')$, the bulk propagator can be decomposed into two distinct terms

$$ G_{++}^{\lambda}(s; \tau', \tau'') = G_{-+}^{\lambda}(s; \tau', \tau'') + [G_{+-}^{\lambda}(s; \tau'', \tau') - G_{-+}^{\lambda}(s; \tau', \tau'')] \Theta(\tau'' - \tau'), \quad (2.31) $$

where the first term is factorised in time, while the second term is still nested. Consequently, the integral $F_{++}^{(\lambda)}$ can be expressed as the sum of two components,
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The full boost-breaking helical seed correlator $F^{(\lambda)}(u, v)$ solved using the bootstrap method, as function of the dimensionless kinematic variable $u$ with fixed $v = 5$, for $\mu = 1$ (left panel) for $\kappa = 0$ (solid line) and $\kappa = 0.5$ (dashed line), and $\mu = 3$ (right panel) for $\kappa = 0$ (solid line) and $\kappa = 1$ (dashed line). We always consider the enhanced helicity $\lambda = -1$ when $\kappa \neq 0$. We include terms in the series from $n = 1, \dots, N$. The shaded gray region corresponds to kinematic configurations outside the unit circle which covers the case of a reduced sound speed (here $c_s = 0.1$) for the external field.</figcaption></figure>

$$F_{++}^{(\lambda)} = F_{F}^{(\lambda)} + F_{N}^{(\lambda)}, \text{ that can be schematically represented as}$$
$$= \text{[factorised part]} + \text{[non-factorised part]}. \qquad (2.32)$$

Here, the first term on the right-hand side of (2.31) is the *factorised* part, which is straightforward to evaluate due to the time-integral factorisation property. The second term is *nested*. The remarkable aspect of the bulk cutting rule is that the purest signals are often the simplest. As a direct consequence of microcausality [80, 81], these oscillatory cosmological collider signals, which include both local and non-local types, entirely stem from the factorised piece, while the nested term only contributes to the background, which is analytic in the momentum ratio.

This cutting rule works perfectly in the case of a unit sound speed, even when the chemical potential effects are taken into account. However, it does not apply to certain kinematical regions when we have a small sound speed. To illustrate this more concretely, let us focus on the integral associated with the factorised part, which contains terms of the form

$$F_{F}^{(\lambda)} \supset \int_{-\infty^{+}}^{0} d\tau e^{ic_{s}k_{12}\tau} \sigma_{s}^{\lambda}(\tau) \sim \int_{-\infty^{+}}^{0} d\tau e^{i(c_{s}k_{12}-s)\tau}, \qquad (2.33)$$

where we have evaluated the UV behaviour of the massive mode function $\tau \to -\infty$. When $c_s = 1$, although the integral contains a spurious folded pole, which should be cancelled by the integral $F_N^{(\lambda)}$, the overall integral remains convergent in the physical kinematic regions due to the triangle inequality $k_{12} > s$. However, a small sound
speed changes the situation significantly, as in certain kinematic regions, $c_s k_{12}$ can become smaller than $s$, leading to a divergence when $\tau \to -\infty^+$. What we are doing here is actually the analytic continuation of the bulk cutting rule. Since the term $F_F^{(\lambda)}$ is factorised, it must satisfy the homogeneous boundary equations and contains all the information about the cosmological collider signals. In contrast, $F_N^{(\lambda)}$ satisfies the inhomogeneous part, as the Heaviside step function generates the source term in the bootstrap differential equation. For regions where $c_s k_{12} > s$, (2.24) and (2.32) are in one-to-one correspondence, as all terms are convergent. When $c_s k_{12} < s$, the resummation of the series allows us to analytically continue each term in (2.24) for a small sound speed, ensuring convergence. This is why we should choose the branch $\tilde{u} = \tilde{u} - i\epsilon$ (or equivalently, $s = s - i\epsilon$, $k_{12,34} = k_{12,34} - i\epsilon$) in the analytical expression when $\tilde{u} > 1$. From the bulk time integral, this corresponds to deforming the contour such that in the UV it behaves as $\sim e^{(s-c_s k_{12})\tau}$ when $\tau \to -\infty$, then it is convergent for cases with a small sound speed. In other words, this process generalises the bulk cutting rule by analytically continuing each term in the bulk time integrals, which can be represented schematically as:

$$ \mathcal{A} \cdot \mathcal{H} \cdot \mathcal{H} = \text{P.V.} \left( \begin{array}{c} \text{F} \\ \text{×} \\ \text{G} \end{array} \right), \quad \mathcal{S} = \text{P.V.} \left( \begin{array}{c} \text{N} \\ \text{×} \\ \text{G} \end{array} \right), \quad (2.34) $$

where the operation P.V. represents extracting the finite part (Cauchy principal value) of the corresponding bulk integral.

* Parity-odd correlator factorisation: exchanging massive spinning fields with a chemical potential can give rise to a large parity-odd four-point correlator [27]. Remarkably, given the constraints from unitarity and scale invariance, parity-odd parts of correlators must factorise, and the total-energy singularity ($k_T \equiv \sum_i k_i \to 0$) should be absent at tree-level [77, 82, 83]. This result was derived from the properties of bulk-bulk propagators under a set of mild assumptions, which are applicable to spinning fields with a chemical potential and sound speed. With the complete analytical expressions now available, we are in a position to perform an explicit consistency check. Consider the example previously discussed, with the interaction (2.11), the parity-odd component is proportional to

$$ \langle \varphi_{k_1} \varphi_{k_2} \varphi_{k_3} \varphi_{k_4} \rangle'_{\text{PO}} \propto i \hat{s} \cdot (\hat{k}_1 \times \hat{k}_3) (1 - k_1 \partial_{k_1}) (1 - k_3 \partial_{k_3}) \sum_{\substack{a,b=\pm \\ \lambda=\pm 1}} \lambda [\mathcal{D}_{k_{12}}^a]^3 [\mathcal{D}_{k_{34}}^b]^3 F_{ab}^{(\lambda)}. \quad (2.35) $$

Given that the parity-odd part is purely imaginary, we have omitted certain real prefactors for the sake of simplicity. The homogeneous solutions are already factorised and do not exhibit any $k_T$-singularity. The theorem thus asserts that the contribution from the particular solution to the parity-odd correlator must vanish, as it originates from the fully connected part of the bulk integral, and also contains the $k_T$-singularity. To verify this, we note that the particular solution (2.29) satisfies the following property

$$ S_{+}^{(\lambda)} = [S_{-}^{(\lambda)}]^* = S_{-}^{(-\lambda)}. \quad (2.36) $$
Then, the summation over helicity gives

$$ \sum_{\substack{a=\pm \\ \lambda=\pm 1}} \lambda S_a^{(\lambda)} = S_+^{(+1)} + S_-^{(+1)} - S_+^{(-1)} - S_-^{(-1)} = 0, \qquad (2.37) $$

which is indeed consistent with the parity-odd factorisation theorem.

* *Comparison with the cosmological low-speed collider result:* as mentioned above, there is an alternative approach to obtaining the convergent solution besides the resummation method employed here. For example, in previous studies involving reduced sound speed but not chemical potential, the ansatz for the particular solution was formulated as $\sum_{m,n=0}^{\infty} (a_{mn} + b_{mn} \log(u)) u^{-m} (u/v)^n$ when $u < v$ [43]. The negative power of $u$ here manifestly ensures the convergence of series when $v > u > 1$. We can adopt a modified ansatz similar to those used in previous studies. However, deriving and solving the recurrence relations for the coefficients proves to be non-trivial. This raises the question: how are these two different methods connected? To facilitate a comparison between both methods, we set the chemical potential to zero, $\kappa = 0$. In this case, the solutions for the particular case with unit sound speed reads

$$ S_+ = \sum_{m,n=0}^{\infty} \frac{(-1)^n (n+1)_{2m}}{2^{2m+1} \left(\frac{n}{2} + \frac{1}{4} + \frac{i\mu}{2}\right)_{m+1} \left(\frac{n}{2} + \frac{1}{4} - \frac{i\mu}{2}\right)_{m+1}} u^{2m+1} \left(\frac{u}{v}\right)^n. \qquad (2.38) $$

This series is convergent for $u < 1$, but we can still perform the resummation along the $m$-direction. After resummation, the series becomes

$$ S_+ = 2u \sum_{n=0}^{\infty} \frac{(-1)^n}{(n + \frac{1}{2})^2 + \mu^2} \left(\frac{u}{v}\right)^n {}_3F_2 \left[ \frac{1}{4} + \frac{n+1}{2}, \frac{n}{2} + 1; \frac{n}{2} + \frac{i\mu}{2}; u^2 \right]. \qquad (2.39) $$

The two layers of infinite series are now reorganised into a single series involving a hypergeometric function. It can be explicitly checked that the result remains valid even for $c_s < 1$. To connect to the results derived in [43], one can expand the hypergeometric function around $u \to \infty$. This reveals a structure consistent with $[a + b \log(u)]u^{-m}$, where the coefficients involve harmonic numbers. In summary, we started here with a series ansatz expansion around $u \to 0$ but after summation, the result converges at the significantly larger radius. The previously obtained result corresponds to the $u \to \infty$ expansion of this expression. The correspondence is quite natural, as the ansatz $[a + b \log(u)]u^{-m}$ was originally formulated based on an analysis of the bootstrap equations in the limit $u, v \to \infty$.

## 2.3 Via spectral representation

We now turn to the computation of the helical seed correlator using the spectral approach recently developed in [56]. Essentially, we extend the spectral representation of the time-ordered bulk-bulk SK propagator to include a helical chemical potential. This allows us to factorise the time integrals at the expense of introducing an additional spectral integration over the mass parameter of the exchanged field. Performing the remaining spectral integral naturally yields a single series that converges in all physical kinematic configurations, in
particular outside the unit circle for $c_s < 1$. Unlike the bootstrap approach, which relies on a
fortunate mathematical trick to resum a series layer and analytically continue the solution
beyond its convergence disk, the spectral method achieves this continuation intrinsically.

**Spectral representation of helical massive propagators.** The essence of the spectral approach lies in the ability to trade time ordering for a contour integral in the complex plane. Instead of implementing the spectral approach for the time-ordered SK propagator $G^\lambda_{++}$, we introduce the following modified bulk-to-bulk propagator$^5$

$$
\mathcal{G}^{\lambda}(k; \tau', \tau'') = \mathcal{G}_{++}^{\lambda}(k; \tau', \tau'') - \Delta\mathcal{G}^{\lambda}(k; \tau', \tau''), \qquad (2.40)
$$

where

$$
\Delta \mathcal{G}^{\lambda}(k; \tau', \tau'') \equiv \frac{\sigma_k^{\lambda}(\tau_0)}{\sigma_k^{\lambda*}(\tau_0)} \sigma_k^{\lambda*}(\tau') \sigma_k^{\lambda*}(\tau''), \qquad (2.41)
$$

with $\tau_0 \to 0$ being a small late-time cutoff. The additional piece $\Delta\mathcal{G}^\lambda$, symmetric under $\tau' \leftrightarrow \tau''$, enforces the correct boundary conditions as it makes $\mathcal{G}^\lambda$ vanish when either of its time arguments is taken to zero. Also note that it is already factorised in time, so that the resulting time integrals, once the propagator is inserted inside a correlator, are performed easily. After Wick rotating time by defining $\tau \equiv ie^{i\epsilon}\chi$ with $\epsilon > 0$, the integral representation of the bulk-bulk propagator $\mathcal{G}^\lambda$ for a given helicity reads

$$
\mathcal{G}^{\lambda}(k; \chi', \chi'') = \frac{1}{2k} \int_{-\infty}^{+\infty} d\rho \mathcal{N}_{\rho, \kappa}^{\lambda} \frac{W_{-i\lambda\kappa, i\rho}(2k\chi')W_{-i\lambda\kappa, i\rho}(2k\chi'')}{\rho^2 + \nu^2}, \quad (2.42)
$$

where we have defined

$$
\mathcal{N}_{\rho, \kappa}^{\lambda} \equiv \frac{\rho}{2\pi^2} \sinh(2\pi\rho)\Gamma(\frac{1}{2} + i\lambda\kappa - i\rho)\Gamma(\frac{1}{2} + i\lambda\kappa + i\rho). \quad (2.43)
$$

Here, to simplify technicalities, we shall temporarily consider the spectral representation for
light fields lying in the complementary series representation of the de Sitter group, where the
dimensionless mass parameter is given by $\nu = \sqrt{9/4 - m^2/H^2} \in \mathbb{R}$. Since the correlators are
analytic functions of the mass parameter at $\nu = 0$, we are allowed to analytically continue
to heavy masses via $\nu \to i\mu$ at the end of the spectral calculation. In the limit $\kappa \to 0$, we
recover the de Sitter density of states $\mathcal{N}_{\rho,\kappa \to 0}^{\lambda} = \frac{\rho}{\pi} \sinh(\pi\rho)$. The appearance of an imaginary
component when $\kappa \neq 0$ reflects parity violation. The detailed derivation of eq. (2.42) is given
in the insert below. In the original spectral representation of the Feynman propagator for
massive fields (with $\kappa = 0$) [56], the projection of incoming (Hankel function) states onto
outgoing (Bessel) states — which encodes particle production — is effectively performed
with a suitable $i\epsilon$ prescription. Interestingly, here, this projection is captured by combining
the Feynman propagator $G_{++}^{\lambda}$ with the boundary term $\Delta G^{\lambda}$ and considering off-shell heavy
states that are set on shell on the imaginary axis, yielding a propagator of light fields.

<footnote>$^5$This propagator is the one used in computing wavefunction coefficients, see e.g. [77, 84, 85] for more details. For fields in the principal series, the term $\sigma_k^\lambda(\tau_0)/\sigma_k^{\lambda*}(\tau_0)$ oscillates as a function of $\tau_0$. We bypass this complication by first considering fields in the complementary series, for which $\sigma_k^\lambda(\tau_0)/\sigma_k^{\lambda*}(\tau_0)$ becomes constant, and then analytically continue the results.</footnote>
Derivation. In this insert, we prove the spectral representation of the bulk-bulk propagator (2.42), following [86]. We first consider the case $\chi' < \chi''$. Using the connection formula

$$W_{-i\lambda\kappa,i\rho}(z) = \frac{\Gamma(-2i\rho)}{\Gamma(\frac{1}{2}-i\rho+i\lambda\kappa)} M_{-i\lambda\kappa,i\rho}(z) + \frac{\Gamma(2i\rho)}{\Gamma(\frac{1}{2}+i\rho+i\lambda\kappa)} M_{-i\lambda\kappa,-i\rho}(z), \quad (2.44)$$

to project the term $W_{-i\lambda\kappa,i\rho}(2k\chi')$ onto a linear combination of Whittaker-M functions, along with the symmetry relation $W_{-i\lambda\kappa,i\rho}(z) = W_{-i\lambda\kappa,-i\rho}(z)$, the integral (2.42) can be written

$$\frac{-i}{2k} \int_{-\infty}^{+\infty} \frac{d\rho}{\pi} \frac{\rho}{\rho^2 + \nu^2} \frac{\Gamma(\frac{1}{2} + i\lambda\kappa - i\rho)}{\Gamma(1 - 2i\rho)} M_{-i\lambda\kappa, -i\rho}(2k\chi') W_{-i\lambda\kappa, -i\rho}(2k\chi''). \quad (2.45)$$

We have also used the Euler's reflection formula $\Gamma(1 - 2i\rho)\Gamma(2i\rho) = -i\pi/\sinh(2\pi\rho)$ to arrive at this expression. We choose to close the integration contour in the upper-half complex plane, since the arc at infinity vanishes for $\chi' < \chi''$. The integrand exhibits two distinct types of poles: (i) on-shell poles located at $\rho = \pm i\nu$ coming from the term $1/(\rho^2 + \nu^2)$, and (ii) a tower of helical poles located at $\rho = -i(n + \frac{1}{2}) + \lambda\kappa$ coming from the term $\Gamma(\frac{1}{2} + i\lambda\kappa - i\rho)$, that are shifted away from the imaginary axis by the non-zero chemical potential. The analytic structure of the integrand in the complex $\rho$ plane is shown in the following figure

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

By Cauchy's residue theorem, since we close the contour upwards for $\chi' < \chi''$, we pick up only the pole $\rho = +i\nu$, which yields

$$\mathcal{G}^{\lambda}(k; \chi', \chi'') = \frac{1}{2k} \frac{\Gamma(\frac{1}{2} + \nu + i\lambda\bar{\kappa})}{\Gamma(1 + 2\nu)} M_{-i\lambda\kappa, \nu}(2k\chi') W_{-i\lambda\kappa, \nu}(2k\chi''). \quad (2.46)$$

This is not yet the desired expression. To connect with the bulk-bulk propagator $\mathcal{G}^{\lambda}$, we note that

$$\Delta \mathcal{G}^{\lambda}(k; \tau', \tau'') = -e^{i\pi(\nu + \frac{1}{2})} \frac{\Gamma(\frac{1}{2} + \nu + i\lambda\kappa)}{\Gamma(\frac{1}{2} + \nu - i\lambda\kappa)} \sigma_k^{\lambda*}(\tau') \sigma_k^{\lambda*}(\tau''). \quad (2.47)$$

This formula is strictly valid only for light fields as we have dropped the most decaying mode when approaching $\tau_0 \to 0$. For heavy fields $\mu = -i\nu > 0$, both modes $(k\tau)^{\frac{1}{2}\pm i\mu}$ are decaying with the same rate while oscillating. We will see however that the derived final expression for the helical seed correlator can straightforwardly be analytically continued to the heavy
case. In deriving (2.47), we have used $W_{i\lambda\kappa,\nu}^*(z) = W_{-i\lambda\kappa,\nu}(z^*)$, valid for all $z \in \mathbb{C}\setminus\{\mathbb{R}_-\}$ away from the branch cut. By using the (inverse) connection formula

$$ \frac{1}{\Gamma(1+2\nu)} M_{i\lambda\kappa,\nu}(z) = \frac{e^{\pm i\pi(i\lambda\kappa - \nu - \frac{1}{2})}}{\Gamma(\frac{1}{2} + \nu + i\lambda\kappa)} W_{i\lambda\kappa,\nu}(z) + \frac{e^{\mp i\pi\lambda\kappa}}{\Gamma(\frac{1}{2} + \nu - i\lambda\kappa)} W_{-i\lambda\kappa,\nu}(e^{\pm i\pi z}), \quad (2.48) $$

the bulk-bulk propagator for a given helicity reduces to

$$ \mathcal{G}^{\lambda}(k; \tau', \tau'') = \frac{1}{2k} \frac{\Gamma(\frac{1}{2} + \nu + i\lambda\tilde{\kappa})}{\Gamma(1 + 2\nu)} M_{-i\lambda\kappa,\nu}(-2ik\tau') W_{-i\lambda\kappa,\nu}(-2ik\tau'') \Theta(\tau' - \tau'') + (\tau' \leftrightarrow \tau''), \quad (2.49) $$

where we have used the formula $M_{-i\lambda\kappa,\nu}(e^{\pm i\pi z}) = e^{\pm i\pi(\nu + \frac{1}{2})} M_{i\lambda\kappa,\nu}(z)$ to rotate the argument without crossing the branch cut which lies on the negative real axis. By defining $\tau \equiv ie^{i\epsilon}\chi$ with $\epsilon > 0$, the propagator $\mathcal{G}^\lambda$ reads

$$ \mathcal{G}^{\lambda}(k; \chi', \chi'') = \frac{1}{2k} \frac{\Gamma(\frac{1}{2} + \nu + i\lambda\tilde{\kappa})}{\Gamma(1 + 2\nu)} M_{-i\lambda\kappa,\nu}(2k\chi') W_{-i\lambda\kappa,\nu}(2k\chi'') \Theta(\chi'' - \chi') + (\chi' \leftrightarrow \chi''), \quad (2.50) $$

in the Wick rotated time domain, so that the arguments of the Whittaker functions lie on the real axis. This exactly boils down to the form (2.46) for $\chi' < \chi''$. For the case $\chi'' < \chi'$, one needs to expand the second Whittaker term $W_{-i\lambda\kappa,ip}(2k\chi'')$ in (2.42) and close the contour in the lower-half complex plane. This completes the proof.

**Time integrals.** We have now the necessary ingredient to compute the helical seed correlator (2.9). The non-time-ordered contribution $F_{+-}^{(\lambda)}$, for which the time integrals factorise, is given by<fnref n="6" />

$$ F_{+-}^{(\lambda)} = \frac{\pi^2 e^{-\pi \lambda \kappa}}{2 \cosh^2(\pi \mu)} {}_2\tilde{F}_1 \left[ \begin{matrix} \frac{1}{2} - i\mu & \frac{1}{2} + i\mu \\ 1 + i\lambda\kappa & \frac{1}{2} - \frac{1}{2u} \end{matrix} \right] {}_2\tilde{F}_1 \left[ \begin{matrix} \frac{1}{2} - i\mu & \frac{1}{2} + i\mu \\ 1 - i\lambda\kappa & \frac{1}{2} - \frac{1}{2v} \end{matrix} \right], \quad (2.51) $$

where we have used the kinematic variables defined in (2.17). Naturally, we have $F_{+-}^{(\lambda)} = F_{+-}^{(\lambda)*}$. After Wick rotating time, the non-time ordered contribution reads

$$ \begin{aligned} F_{++}^{(\lambda)} &= -s \int_0^{+\infty} \frac{d\chi'}{\chi'} \frac{d\chi''}{\chi''} e^{-c_s k_{12} \chi'} \mathcal{G}_{++}^{\lambda}(s; \chi', \chi'') e^{-c_s k_{34} \chi''} \\ &\equiv F_{\mathcal{G}}^{(\lambda)} + F_{\Delta\mathcal{G}}^{(\lambda)}, \end{aligned} \quad (2.52) $$

where we have split the integral according to (2.40). The piece $F_{\Delta\mathcal{G}}^{(\lambda)}$ coming from the difference between the Feynman propagator and the wavefunction bulk-bulk propagator is factorised in time and can be easily computed, yielding

$$ F_{\Delta\mathcal{G}}^{(\lambda)} = \frac{i\pi^2 e^{-\pi(\mu+\lambda\kappa)}}{2\cosh^2(\pi\mu)} \frac{\Gamma(\frac{1}{2}+i\lambda\kappa+i\mu)}{\Gamma(\frac{1}{2}-i\lambda\kappa+i\mu)} {}_2\tilde{F}_1 \left[ \begin{matrix} \frac{1}{2}-i\mu & \frac{1}{2}+i\mu \\ 1+i\lambda\kappa & \frac{1}{2} \end{matrix} \right] {}_2\tilde{F}_1 \left[ \begin{matrix} \frac{1}{2}-i\mu & \frac{1}{2}+i\mu \\ 1+i\lambda\kappa & \frac{1}{2} \end{matrix} \right]. \quad (2.53) $$

<footnote n="6">⁶Here we have analytically continued back to heavy masses via $\nu \to i\mu$. We are guaranteed to land on a correct result since the correlators are analytic functions of the dimensionless mass parameter near the origin $\nu = 0$.</footnote>
Using the spectral representation (2.42) and after swapping the spectral and the time integrals, the remaining contribution reads

$$F_G^{(\lambda)} = \frac{-\pi^2}{2} \int_{-\infty}^{+\infty} \frac{d\rho \mathcal{N}_{\rho,\kappa}^{\lambda}}{\cosh^2(\pi\rho)} \frac{2\tilde{F}_1\left[\frac{1}{2}-i\rho, \frac{1}{2}+i\rho; \frac{1}{2}-\frac{1}{2u}\right] 2\tilde{F}_1\left[\frac{1}{2}-i\rho, \frac{1}{2}+i\rho; \frac{1}{2}-\frac{1}{2v}\right]}{\rho^2 + \nu^2} . \quad (2.54)$$

This spectral representation makes it manifest that the integral $F_G^{(\lambda)}$ is obtained by gluing two off-shell three-point functions, and we naturally expect the solution to be of higher transcendentality than the factorised integral $F_{+-}^{(\lambda)}$ due to the additional spectral integration.

**Collecting poles.** Computing the helical seed correlator boils down to evaluating the spectral integral (2.54). The off-shell three-point function can be written in a compact form using the Legendre function<fnref n="7" />

$$P_{i\rho-1/2}^{-i\lambda\kappa}(z) = \left(\frac{z+1}{z-1}\right)^{-i\lambda\kappa/2} {}_2F_1\left[\frac{1}{2}-i\rho, \frac{1}{2}+i\rho; \frac{1}{2}; \frac{z}{2}\right], \quad (2.55)$$

so that $F_G^{(\lambda)}$ can be written as

$$\begin{aligned} F_G^{(\lambda)} &= \frac{-\pi^2}{2} \left[ \left( \frac{1-u}{1+u} \right) \left( \frac{1-v}{1+v} \right) \right]^{-\lambda i \kappa / 2} \int_{-\infty}^{+\infty} \frac{d\rho \mathcal{N}_{\rho,\kappa}^{\lambda}}{\cosh^2(\pi\rho)} \frac{P_{i\rho-1/2}^{-i\lambda\kappa}(u^{-1}) P_{i\rho-1/2}^{-i\lambda\kappa}(v^{-1})}{\rho^2 + \nu^2} \\ &\equiv \left[ \left( \frac{1-u}{1+u} \right) \left( \frac{1-v}{1+v} \right) \right]^{-\lambda i \kappa / 2} \tilde{F}_G^{(\lambda)}, \end{aligned} \quad (2.56)$$

where we have defined the “stripped” integral $\tilde{F}_G^{(\lambda)}$. We finish this last integral using Cauchy residue’s theorem. To this end, we need to understand the analytic structure of the meromorphic integrand, i.e. identify poles, and disentangle the large-$\rho$ behaviour to close the contour accordingly. The Legendre-$P$ basis behaves well (meaning that it has a single mode) around the origin $\rho \approx 0$, whereas the associated Legendre-$Q$ basis behaves well at infinity $\rho \to \infty$. We therefore use the connection formula

$$\Gamma(\frac{1}{2} + i\rho + \lambda i\kappa)\Gamma(\frac{1}{2} - i\rho + \lambda i\kappa)P_{i\rho-1/2}^{-\lambda i\kappa}(x) = \frac{e^{-i\pi/2 + \pi\lambda\kappa}}{\sinh(\pi\rho)} \left[Q_{-i\rho-1/2}^{\lambda i\kappa}(x) - Q_{+i\rho-1/2}^{\lambda i\kappa}(x)\right], \quad (2.57)$$

to decompose the integral into two pieces $\tilde{F}_G^{(\lambda)} = \tilde{F}_G^{(\lambda),0} + \tilde{F}_G^{(\lambda),<,>}$, with (we have accounted for the shadow symmetry of the integrand $\mu \leftrightarrow -\mu$)

$$\begin{aligned} \tilde{F}_G^{(\lambda),0} &\equiv \frac{2e^{2\pi\lambda\kappa}}{\pi} \int_{-\infty}^{+\infty} \frac{d\rho}{2\pi} \frac{|\Gamma(1+i\rho)|^2 |\Gamma(\frac{1}{2}+i\rho)|^2}{\Gamma(\frac{1}{2}+i\lambda\kappa+i\rho)\Gamma(\frac{1}{2}+i\lambda\kappa-i\rho)} \frac{Q_{-i\rho-1/2}^{i\lambda\kappa}(u^{-1})Q_{-i\rho-1/2}^{i\lambda\kappa}(v^{-1})}{\rho^2 + \nu^2}, \\ \tilde{F}_G^{(\lambda),<,>} &\equiv -\frac{2e^{2\pi\lambda\kappa}}{\pi} \int_{-\infty}^{+\infty} \frac{d\rho}{2\pi} \frac{|\Gamma(1+i\rho)|^2 |\Gamma(\frac{1}{2}+i\rho)|^2}{\Gamma(\frac{1}{2}+i\lambda\kappa+i\rho)\Gamma(\frac{1}{2}+i\lambda\kappa-i\rho)} \frac{Q_{-i\rho-1/2}^{i\lambda\kappa}(u^{-1})Q_{+i\rho-1/2}^{i\lambda\kappa}(v^{-1})}{\rho^2 + \nu^2}, \end{aligned} \quad (2.58)$$

<footnote n="7">⁷Note that this function is implemented in Mathematica by LegendreP[iμ - 1/2, -iλκ, 3, z] where the “3” encodes the correct branch cut discontinuity.</footnote>
where we have used the Legendre duplication formula $\Gamma(z)\Gamma(\frac{1}{2} + z) = 2^{1-2z}\sqrt{\pi}\Gamma(2z)$.<fnref n="8" />
Assuming $|u| < |v|$, we close the contour in the upper-half complex plane for both integrals as again the arcs at infinity vanish. The two integrands have several poles with different physical origins:

* The particle production poles coming from $1/(\rho^2 + \nu^2)$ are located at $\rho = \pm i\nu$. Picking up these poles extracts the non-analytic signal contributions, which read

$$ 
\begin{align} 
\tilde{F}_{G}^{(\lambda),0} & \supset -\frac{ie^{2\pi\lambda\kappa}}{\pi\mu} \frac{|\Gamma(\frac{1}{2} + i\mu)|^2 |\Gamma(1 + i\mu)|^2}{\Gamma(\frac{1}{2} + i\lambda\kappa + i\mu)\Gamma(\frac{1}{2} + i\lambda\kappa - i\mu)} Q_{i\mu-1/2}^{i\lambda\kappa}(u^{-1}) Q_{i\mu-1/2}^{i\lambda\kappa}(v^{-1}), \\ 
\tilde{F}_{G}^{(\lambda),<,>} & \supset \frac{ie^{2\pi\lambda\kappa}}{\pi\mu} \frac{|\Gamma(\frac{1}{2} + i\mu)|^2 |\Gamma(1 + i\mu)|^2}{\Gamma(\frac{1}{2} + i\lambda\kappa + i\mu)\Gamma(\frac{1}{2} + i\lambda\kappa - i\mu)} Q_{i\mu-1/2}^{i\lambda\kappa}(u^{-1}) Q_{-i\mu-1/2}^{i\lambda\kappa}(v^{-1}), 
\end{align} 
\tag{2.59} $$

where we have set $\nu = i\mu$. These expressions are factorised in $u$ and $v$ as they should.

* The tower of density of states poles coming from the factor $|\Gamma(1 + i\rho)|^2$ are located on the entire imaginary axis $\rho = \pm i(n + 1)$. These two towers exactly cancel again each other, as also noted in [56].

* The EFT poles coming from the factor $|\Gamma(\frac{1}{2} + i\rho)|^2$ are located at $\rho = \pm i(n + \frac{1}{2})$. We therefore pick up the poles in the upper-half complex plane and sum over $n \ge 0$. Writing the Legendre-Q function in its hypergeometric form<fnref n="9" />

$$ 
\begin{align} 
Q_{-i\rho-1/2}^{i\lambda\kappa}(z) & \equiv e^{-\pi\lambda\kappa} \sqrt{\frac{\pi}{2}} \Gamma(\frac{1}{2} - i\rho + \lambda i\kappa) (z^2 - 1)^{\lambda i\kappa/2} 2^{i\rho} z^{-\frac{1}{2} + i\rho - \lambda i\kappa} \nonumber \\ 
& \quad \times {}_2F_1 \left[ \frac{3/2 - i\rho + \lambda i\kappa}{2} \quad \frac{1/2 - i\rho + \lambda i\kappa}{2} \quad \frac{1}{z^2} \right], 
\end{align} 
\tag{2.60} $$

the EFT contribution coming from the symmetric integral $F_G^{(\lambda),0}$ reads

$$ 
\begin{align} 
F_{G, \text{EFT}}^{(\lambda),0} = & -\pi \left[(u+1)(v+1)\right]^{\lambda i \kappa} \times \sum_{n=0}^{+\infty} \frac{(1+2n)}{(n+\frac{1}{2})^2 + \mu^2} \left(\frac{v}{2}\right)^{n+1} \left(\frac{u}{2}\right)^{n+1} \nonumber \\ 
& \times \frac{\Gamma(1+n+i\lambda\kappa)}{\Gamma(-n+i\lambda\kappa)} {}_2F_1 \left[ \frac{1+n}{2} + \frac{\lambda i \kappa}{2} \quad 1+\frac{n}{2} + \frac{\lambda i \kappa}{2} ; u^2 \right] {}_2F_1 \left[ \frac{1+n}{2} + \frac{\lambda i \kappa}{2} \quad 1+\frac{n}{2} + \frac{\lambda i \kappa}{2} ; v^2 \right]. 
\end{align} 
\tag{2.61} $$

This EFT series vanishes for $\kappa = 0$ due to the factor $1/\Gamma(-n + i\lambda\kappa)$, in accordance with [56], meaning that this EFT series is a purely boost-breaking contribution to the full correlator. The second EFT series coming from $F_G^{(\lambda),<,>}$ reads

$$ 
\begin{align} 
F_{G, \text{EFT}}^{(\lambda),<,>} = & u \left[(u+1)(v+1)\right]^{\lambda i \kappa} \times \sum_{n=0}^{+\infty} \frac{(-1)^n}{(n+\frac{1}{2})^2 + \mu^2} \left(\frac{u}{v}\right)^n \nonumber \\ 
& \times {}_2F_1 \left[ \frac{1+n}{2} + \frac{\lambda i \kappa}{2} \quad 1+\frac{n}{2} + \frac{\lambda i \kappa}{2} ; u^2 \right] {}_2F_1 \left[ \frac{1-n}{2} + \frac{\lambda i \kappa}{2} \quad -\frac{n}{2} + \frac{\lambda i \kappa}{2} ; v^2 \right]. 
\end{align} 
\tag{2.62} $$

<footnote n="8">⁸Note that whenever we write $|\Gamma(x + iy)|^2$ with both real $x$ and $y$, we actually mean $\Gamma(x + iy)\Gamma(x - iy)$ as there is an ambiguity when locating the poles.</footnote>

<footnote n="9">⁹Note that the factor $\Gamma(\frac{1}{2} - i\rho + \lambda i\kappa)$ cancels against itself inside the integral. As such, this term does not bring additional tower of poles.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The full boost-breaking helical seed correlator $F^{(\lambda)}(u, v)$ solved using the spectral representation method, as function of the dimensionless kinematic variable $u$ with fixed $v = 5$, for $\mu = 1$ (left panel) for $\kappa = 0$ (solid line) and $\kappa = 0.5$ (dashed line), and $\mu = 3$ (right panel) for $\kappa = 0$ (solid line) and $\kappa = 1$ (dashed line). We by default consider the enhanced helicity $\lambda = -1$ when $\kappa \neq 0$. We include terms in the series from $n = 1, \dots, N$. The shaded gray region corresponds to kinematic configurations outside the unit circle which covers the case of a reduced sound speed (here $c_s = 0.1$) for the external field.</figcaption></figure>

This expression reduces to the one found in [56] for $\kappa = 0$. The total EFT series contribution is given by the sum of these two pieces: $F_{G,EFT}^{(\lambda)} = F_{G,EFT}^{(\lambda),0} + F_{G,EFT}^{(\lambda),<,>}$. The result for $|u| > |v|$ is simply found after swapping $u$ and $v$.

**Full result.** The full result is found by summing all contributions, namely the factorised piece (2.51), the boundary term (2.53), the collider contributions from the time-ordered integrals (2.59) together with (2.56), as well as the two EFT series (2.61) and (2.62). We have implemented the result in the Mathematica notebook available on the Github repository linked in [93]. We have explicitly checked that this spectral representation exactly matches the expression found using the bootstrap approach, in all kinematic configurations. A noticeable feature of the spectral representation is that the solution is already partially resummed, yielding a single series of (a product of) hypergeometric functions. This automatically encodes the analytic continuation of the helical seed correlator outside the unit kinematic circle. In figure 4, we show the correlator $F^{(\lambda)}$ as function of the kinematic variables extended outside the unit circle up to $u = 10$ (which describes the case $c_s = 0.1$). We notice that the solution converges well and stabilise when the series contains $N \sim O(50)$ terms. However, compared to the bootstrap representation shown in figure 3, convergence is much slower, and the solution when summing just a few terms in the series features oscillating behaviours before converging to the exact result. We also notice a more pronounced discontinuity as the kinematic variable $u$ crosses the unit circle $u = 1$, which fades away as we increase the number of terms in the series.

**Discussion.** Further comments are in order.

* *Boost-breaking EFT series:* the spectral method naturally yields a partially resummed representation. This results in the two EFT series (2.61) and (2.62). Physically, these
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Amplitude of the series coefficients $|F_n^{(\lambda)}|$ to the boost-breaking helical seed correlator $\sum_n F_n^{(\lambda)}$ as function of $n$ for the bootstrap representation (in blue) and the spectral representation (in red), both inside the unit kinematic circle (left panel) and outside (right panel). We have chosen $\mu = 1, \kappa = 0.5$ (upper panel) and $\mu = 5, \kappa = 2$ (lower panel). We consider by default the negative helicity $\lambda = -1$ when $\kappa \neq 0$.</figcaption></figure>

series correspond to summing over the quasi-normal modes of the mediated field. To see this, let us expand the helical mode function (2.7) in powers of $(-k\tau)$ at late times

$$ \sigma_k^\lambda(\tau) = \frac{e^{-\pi\lambda\kappa/2}}{\sqrt{2k}} \sum_{n=0}^{+\infty} \sum_{\pm} c_n^\pm (-k\tau)^{\frac{1}{2}+n\pm i\mu}, \quad (2.63) $$

with

$$ c_n^\pm = \pm \frac{e^{i\frac{\pi}{2}(n+\frac{1}{2})}}{n!} \frac{\pi 2^{\frac{1}{2}\pm i\mu}}{\Gamma(\frac{1}{2}\mp i\mu - i\lambda\kappa)} {}_2F_1 \left[ -n \quad \frac{1}{2} \pm i\mu - i\lambda\kappa \quad ; \ 2 \right] \frac{e^{\pm\pi\mu/2}}{\sinh(2\pi\mu)}. \quad (2.64) $$

This series is the gradient expansion of the mode function as $(-k\tau)^n = k_p^n$, with $k_p \equiv k/a$ being the physical redshifting momentum (we have set $H = 1$). A noticeable feature is that this expansion contains odd powers of $k_p$, which reflects the breaking of boosts by the chemical potential. We can therefore interpret the EFT series $F_{G,EFT}^{(\lambda),0}$ as summing over odd quasi-normal modes, and the EFT series $F_{G,EFT}^{(\lambda),<,>}$ as summing over the even quasi-normal modes. Therefore, when restoring boosts by setting $\kappa = 0$, we naturally expect the entire tower of odd quasi-normal modes to vanish, and consequently $F_{G,EFT}^{(\lambda),0} = 0$ at the level of the helical seed correlator. This is consistent with the fact that when the chemical potential vanishes, there is no linear gradient term in the dispersion relation of the massive field, as can be seen from the equation of motion (2.6). Explicitly, setting $\kappa = 0$ in the coefficients (2.64), we have [87]

$$ {}_2F_1 \left[ -n \quad \frac{1}{2} \pm i\mu \quad ; \ 2 \right] = \frac{\Gamma\left(\frac{n}{2} + \frac{1}{2}\right) \Gamma(1 \pm i\mu)}{2\sqrt{\pi} \Gamma\left(1 + \frac{n}{2} \pm i\mu\right)} \times \left[(-1)^n + 1\right], \quad (2.65) $$

so that all odd powers of $(-k\tau)$ vanish.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** CPU time in milliseconds to evaluate the boost-breaking helical seed correlator $F^{(\lambda)}(u, v)$ including terms in the series from $n = 1, \dots, N$ as function of $N$ for the bootstrap representation (in blue) and the spectral representation (in red), both inside the unit kinematic circle (left panel) and outside (right panel). We have chosen $\mu = 1, \kappa = 0.5$ (upper panel) and $\mu = 5, \kappa = 2$ (lower panel). We consider by default the negative helicity $\lambda = -1$ when $\kappa \neq 0$. These computations were performed with Mathematica on a Macbook Pro with M1 CPU and running MacOS 11.6 Big Sur.</figcaption></figure>

* *Comparison with bootstrap result:* both the spectral and the bootstrap methods yield a convergent series representation albeit with different convergence rates. In figure 5, we show the amplitude of the series coefficients $|F_n^{(\lambda)}|$ of the boost-breaking helical seed correlator $\sum_n F_n^{(\lambda)}$ as function of $n$. Both inside and outside the kinematic unit circle, the bootstrap series has a faster convergence rate. Keeping track of the mass parameter $\mu$ dependence only, this is because each term scales as $F_n^{(\lambda)} \sim 1/\mu^{2n}$, as can be seen from (2.29), whereas the spectral series scales as $F_n^{(\lambda)} \sim 1/(n^2 + \mu^2)$, see (2.62). A noticeable feature of the spectral representation series is that the coefficient amplitude decrease $|F_n^{(\lambda)}|$ is not monotonic as $n$ increases outside the unit kinematic circle, but exhibits an (slowly decreasing) oscillating behaviour. This causes the solution to oscillate in the $u$-domain as, for a fixed number of terms in the series $n$, some coefficients at larger $n$ would be dominant. Interestingly, inside the unit circle, both series' evaluation times are of the same order of magnitude $\sim O(10\text{ms})$, as can be seen in figure 6. However, beyond the unit kinematic circle, the computational cost of the bootstrap representation increases linearly, whereas the spectral representation remains stable. This difference arises because the analytic continuation of the $_2F_1$ hypergeometric function is significantly more efficient than that of $_3F_2$. However, since the bootstrap series converges more rapidly, an accurate result can be achieved with only a few terms, whereas the spectral representation requires summing substantially more terms to reach a comparable level of accuracy.
## 3 Classifying non-analyticities

We now introduce the saddle point method supplemented with the WKB approximation, as a powerful tool to classify and gain physical insight into non-analytic signals of cosmological correlators. This approach relies solely on analysing the location of saddle points in the complex time plane — dominating the time integrals—, providing a natural explanation for the behaviour of these signals by tracking the motion of these saddles driven by external kinematics. By evaluating the time integrals using the saddle point method, we obtain approximate solutions for a range of correlators, expressed entirely in terms of elementary functions.

In this section, we start by outlining the general strategy to classify non-analyticities of correlators and approximating cosmological collider signals. Then, we apply the method to concrete cases with increasing complexity, and provide new cosmological collider templates for future cosmological surveys.

### 3.1 Extracting & approximating non-analyticities

Without loss of generality, we work with the exchange four-point seed correlators $F$ defined in section 2.1, whose kinematics is described by external energies, $k_{12}$ and $k_{34}$, and the internal momentum $s$ in the $s$-channel. We include a reduced speed of sound for external legs and a helical chemical potential for the exchanged massive particle. These boost-breaking effects will be turned on in turn lately. Our approximation strategy consists of four steps.

* Apply bulk cutting rules: non-analytic scaling of cosmological correlators — and therefore cosmological collider signals — come from factorised contributions. We limit ourselves to the kinematic region $k_{12} > k_{34}$ and focus on the factorised contribution [51, 80, 81]

$$F_F = \int_{-\infty^+}^0 d\tau' e^{ic_s k_{12} \tau'} \int_{-\infty^+}^0 d\tau'' e^{ic_s k_{34} \tau''} \mathcal{G}_{-+}(s; \tau', \tau''). \quad (3.1)$$

Here, we have removed powers of conformal time. We stress that this choice simplifies our technical discussion later on, but not at the cost of generality nor physicality. Note that we have also dropped overall kinematic factors that make the integral dimensionless. These prefactors can easily be restored at the end. Factorising the Wightman function $\mathcal{G}_{-+} \equiv \mathcal{G}_>^-$ into mode functions, we can separate the integral into two parts $F_F = F_L^{(3)} \times F_R^{(3)}$ with

$$F_L^{(3)} = \int_{-\infty^+}^0 d\tau' e^{ic_s k_{12} \tau'} \tilde{\sigma}_s(\tau'), \quad F_R^{(3)} = \int_{-\infty^+}^0 d\tau'' e^{ic_s k_{34} \tau''} \tilde{\sigma}_s^*(\tau''), \quad (3.2)$$

where $\tilde{\sigma}_s$ is the rescaled mode function of the exchanged field, see eq. (2.8). In most of the cases, these integrals can be easily computed and result in hypergeometric functions.

* Apply the WKB approximation: the rescaled mode function satisfies an equation of motion of the form

$$\tilde{\sigma}_s'' + \omega_s^2(\tau)\tilde{\sigma}_s = 0, \quad (3.3)$$

where $\omega$ is an effective frequency for the $s$-mode that is theory-dependent. In particular, it depends on the momentum magnitude $s$ and other characteristics of the exchange
field such as its mass and the chemical potential. Instead of directly solving the exact mode function, we write the solution as a time-dependent superposition of positive- and negative-frequency parts, both solved via the superadiabatic WKB approximation [68, 88, 89]

$$ \tilde{\sigma}_s(\tau) \approx \alpha_s(\tau)u_s(\tau) + \beta_s(\tau)u_s^*(\tau), \quad \text{with} \quad u_s(\tau) \equiv \frac{e^{-i \int_{\tau_i}^{\tau} d\tau' \omega_s(\tau')}}{\sqrt{2\omega_s(\tau)}}, \qquad (3.4) $$

where $\tau_i$ is an arbitrary reference (conformal) time. The integration contour follows the negative real axis with a slight shift towards the positive imaginary direction due to the $i\epsilon$-prescription. The time-dependent Bogoliubov coefficients are approximately given by

$$ \alpha_s(\tau) \approx 1, \quad \beta_s(\tau) \approx -iS(\tau)e^{-F_s(\tau)}, \qquad (3.5) $$

where

$$ S(\tau) = \frac{1}{2} \left[ 1 + \operatorname{erf} \left( \frac{-\operatorname{Im} F_s(\tau)}{\sqrt{2|\operatorname{Re} F_s(\tau)|}} \right) \right] \qquad (3.6) $$

is the known as the Stokes multiplier and $F_s$ is the Dingle's singulant defined by [88]

$$ F_s(\tau) = -2i \int_{\tau_c}^{\tau} \omega_s(\tau') d\tau', \qquad (3.7) $$

with $\tau_c$ being the turning points such that $\omega_s^2(\tau_c) = 0$ with $\operatorname{Im}\tau_c < 0$. In the early-time limit $\tau \to -\infty$, and the Stokes multiplier asymptotes to zero i.e. $S(-\infty) = 0$, thereby recovering the Bunch-Davies initial condition. However, it asymptotes to unity in the late-time limit, $S(0^-) = 1$, indicating the emergence of the negative-frequency mode and spontaneous particle production. The exact moment of particle production can be read from the Stokes line crossing time $\tau = \tau_*$ where $\operatorname{Im}F_s(\tau_*) = 0$ and $S(\tau_*) = 1/2$ [68, 89].

* Use the saddle point method to find resonant times: plugging the WKB approximation (3.4) in the bulk time integrals $F_L^{(3)}$ and $F_R^{(3)}$, all integrals are of the form

$$ \int_C d\tau g(\tau) e^{\mu f(\tau)}, \qquad (3.8) $$

where $g$ and $f$ are theory-dependent functions of the complex conformal time $\tau$, and $\mu$ denotes the dimensionless mass parameter of the exchanged particle. Here, $C$ denotes the specific contour on the complex plane that encodes the $i\epsilon$ prescription. In the large-mass regime $\mu \gg 1$, integrals of this type are known to be solvable as an asymptotic series in $\mu^{-1}$ via the saddle-point method. Such saddle-point approximations are capable of partially capturing effects exponentially suppressed by the mass $\mu$, and is therefore suitable for understanding non-perturbative physics beyond the local large-mass EFTs. We give details about this method in appendix A. Solving for the saddle point $\tau_\bullet$ such that $f'(\tau_\bullet) = 0$, the integration contour $C$ is continuously deformed to a new contour $C'$ that passes through the saddle point in the steepest direction. Note that in practice, whether or not $C'$ has to pass exactly along the steepest descent direction is usually
unimportant, since the integrand is often analytic around $\tau_•$ and small deformations are allowed. Expanding the integrand around the saddle, the Gaussian integral can be finished, leading to

$$ \int_C d\tau g(\tau) e^{\mu f(\tau)} \approx g(\tau_•) e^{\mu f(\tau_•)} \sqrt{\frac{-2\pi}{\mu f''(\tau_•)}} [1 + \mathcal{O}(\mu^{-1})] . \quad (3.9) $$

Note that higher-order terms are small due to the small variance of the Gaussian. Therefore, at leading order, the resulting integral is fully controlled by the location of the saddle, which is a function of the external kinematics, i.e. $\tau_• = \tau_•(k_{12}, k_{34}, s)$.

* **Classify cosmological collider signals:** to relate the mathematical notion of saddle points to the physical bulk evolution, we note that massive fields are constantly produced and decay during inflation. Due to the lack of time-translational symmetry and hence the non-conservation of energy, particles can be produced either in pairs from the propagator, or in singlet from the interaction vertices, before decaying. These production and decay events happen in the form of resonant transitions that manifest themselves as saddles on the complex time plane. The time spent by the particle between these saddles is thus recorded in the resultant cosmological collider signals. One can also estimate the size of the cosmological signal by inspecting the location of the saddle. The rule of thumb is

Saddle points lying in the physical region i.e. $\tau_• \in \mathbb{R}_- \to \mu^{-\#}$ suppression,

Saddle points lying in the unphysical region i.e. $\tau_• \notin \mathbb{R}_- \to e^{-\#\mu}$ suppression, \hfill (3.10)

where $\#$ denotes an order-one positive number. One can intuitively understand this as an “on-shell” versus “off-shell” effect. To summarise, by studying the kinematic dependence of the saddle point $\tau_• = \tau_•(k_{12}, k_{34}, s)$, we can classify cosmological collider signals into different categories with distinct shapes and sizes.

Through this procedure, we learn a beautiful correspondence between the observational, the mathematical, and the physical aspects of cosmological collider signals: (i) the momentum ratio $k_L/k_S$, where $k_L$ is the momentum of a generic long-wavelength mode and $k_S$ is that of a generic short-wavelength one, reflects the kinematic control (more precisely the internal soft limit) of the boundary cosmological correlator on the observational side, (ii) the location of the complex saddle point $\tau_•$ dictates the dominant sources of the time integrals on the mathematical side, and (iii) the lifetime of the massive particle $\Delta\tau$ traces the origin of the cosmological collider signal on the physical side. Modifying any one of the three variables is equivalent to changing the other two, as represented in the picture below.

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Illustration of the two-sheet Riemann surface for the effective frequency $\omega_s(\tau)$ in eq. (3.11). We have retained the topology of the Riemann surface while removing the pole at the origin $\tau = 0$ for better visualisation. The branches that asymptote to red and cyan correspond to Sheet I and Sheet II, respectively, as defined in (3.12). The curve connecting $\tau_c$ and $\tau_c^*$ is the branch cut, along which the two sheets are smoothly sewn together. However, discontinuity appears when crossing the branch cut and landing on the same sheet. We have cartooned the same “journey of the spider spider” to better demonstrate the structure of the two-sheet Riemann surface.</figcaption></figure>

## 3.2 Application to tree-level exchange correlators

We now apply the strategy outlined above to concrete examples with increasing complexity.

### 3.2.1 de Sitter-invariant scenario

Let us first consider the minimal canonical scenario of an exchanged massive scalar field with zero chemical potential $\kappa = 0$ and with a trivial sound speed $c_s = 1$. The effective frequency in this case is given by

$$ \omega_s^2(\tau) = s^2 + \frac{\mu^2}{\tau^2}, \qquad (3.11) $$

which has a pole at $\tau = 0$ and two simple zeros at $\tau_c, \tau_c^* = \mp i\mu/s$. Taking the square root gives a double-valued function defined on two Riemann sheets $\mathbb{C}_{I,II}$ that we define by

$$ \omega_s(\tau) = +\sqrt{s^2 + \frac{\mu^2}{\tau^2}}, \quad \tau \in \mathbb{C}_I, \quad \text{and} \quad \omega_s(\tau) = -\sqrt{s^2 + \frac{\mu^2}{\tau^2}}, \quad \tau \in \mathbb{C}_II, \qquad (3.12) $$

where $\sqrt{\cdot}$ denotes the standard principal-valued square root function. On either sheets, there is a branch cut that joins $\tau_c, 0$ and $\tau_c^*$, but there are no discontinuities between $\tau_c, 0$ and $\tau_c^*$ across different sheets. Instead, the discontinuities on the imaginary axis are from $\tau_c, \tau_c^*$ to $\mp i\infty$. The corresponding surface is represented in figure 7. We now follow the recipe outlined in the previous section. The first factorised time integral $F_L^{(3)}$ is given by

$$ F_L^{(3)} \approx \int_{-\infty+}^{0} \frac{d\tau'}{\sqrt{2\omega_s(\tau')}} \left[ e^{i(k_{12}\tau' - \int_{\tau_i}^{\tau'} d\tau \omega_s(\tau))} + \beta_s(\tau')e^{i(k_{12}\tau' + \int_{\tau_i}^{\tau'} d\tau \omega_s(\tau))} \right] \qquad (3.13) \\ \equiv F_L^{(3),\alpha} + F_L^{(3),\beta}, $$
where the superscripts α and β keep track of the positive- and negative-frequency parts of the massive mode functions. Here, we recall that $\tau_i$ is an arbitrary reference time which can be taken to be $\tau_i \to -\infty$ without loss of generality. To find the saddle points, we solve

$$ \frac{d}{d\tau'} \left( k_{12}\tau' \mp \int_{\tau_i}^{\tau'} d\tau \omega_s(\tau) \right) = 0. \qquad (3.14) $$

For the integral $F_L^{(3),\alpha}$, the nearest saddle — which contributes the most — lies in the physical domain of $\mathbb{C}_I$

$$ \tau_{\bullet}^{\prime(\alpha)} = \frac{-\mu}{\sqrt{k_{12}^2 - s^2}} \in \mathbb{C}_I. \qquad (3.15) $$

Although $\tau'_\bullet$ and $-\tau'_\bullet$ are both solutions of the saddle-point equation (3.14), according to the rule of thumb, we always pick up the one closest to the physical region since its contribution is the largest. For the integral $F_L^{(3),\beta}$, there are no saddle points in the physical domain. Instead, the nearest saddle point lies in $\mathbb{C}_{II}$

$$ \tau_{\bullet}^{\prime(\beta)} = \frac{+\mu}{\sqrt{k_{12}^2 - s^2}} \in \mathbb{C}_{II}. \qquad (3.16) $$

In addition, the Bogoliubov coefficient is suppressed by $|\beta_s| = e^{-\pi\mu}$. Therefore, we can completely discard this term. In fact, one can show that this term evaluated at the saddle $\tau_{\bullet}^{\prime(\beta)}$ is of order $O(e^{-2\pi\mu})$. Its cousin at $-\tau_{\bullet}^{\prime(\beta)} \in \mathbb{C}_{II}$ is unreachable from the original contour on the first Riemann sheet and therefore does not contribute. For pedagogical reasons, we gather a pedagogical introduction to Riemann surfaces and the saddle point method (as well as concrete worked out examples) in appendix A. Performing the saddle-point approximation for the positive-frequency integral, we obtain

$$ F_L^{(3)} \approx F_L^{(3),\alpha} \approx \frac{\sqrt{\pi\mu}}{(k_{12}^2 - s^2)^{3/4}} e^{-i\pi/4} \exp \left[ -i\mu \operatorname{arccosh} \frac{k_{12}}{s} \right]. \qquad (3.17) $$

We give details of the derivation in appendix A, and represent the complex $\tau'$ plane in figure 8. Note that keeping the dependence on $\tau_i$ gives an unphysical phase that will eventually be cancelled in the final result, after multiplying by $F_R^{(3)}$.

The factorised time integral $F_R^{(3)}$ turns out to be more complicated. Again, let us separate the integral into positive- and negative-frequency modes

$$ F_R^{(3)} \approx \int_{-\infty^+}^0 \frac{d\tau''}{\sqrt{2\omega_s(\tau'')}} \left[ e^{i(k_{34}\tau'' + \int_{\tau_i}^{\tau''} d\tau \omega_s(\tau))} + \beta_s^*(\tau'')e^{i(k_{34}\tau'' - \int_{\tau_i}^{\tau''} d\tau \omega_s(\tau))} \right] \qquad (3.18) \\ \equiv F_R^{(3),\alpha} + F_R^{(3),\beta}. $$

The saddle points now are solutions to the following equation

$$ \frac{d}{d\tau''} \left( k_{34}\tau'' \pm \int_{\tau_i}^{\tau''} d\tau \, \omega_s(\tau) \right) = 0. \qquad (3.19) $$
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Analytic structure of the effective frequency $\omega_s$ in the complex $\tau'$ plane, together with both Riemann sheets $C_I$ and $C_{II}$, for the case $c_s = 1$ and $\kappa = 0$ to compute $F_L^{(3)}$. Here and in subsequent figures, the double lines indicate the disconnection that cannot be crossed between the two sheets. We also represent the turning points $\tau_c$, $\tau_c^*$, and the saddle points $\tau_\bullet^{(\alpha)}$, $\tau_\bullet^{(\beta)}$. The original integration contour is in gray, which encodes the early-time $i\epsilon$ prescription.</figcaption></figure>

The integral $F_R^{(3),\beta}$, despite being suppressed by $|\beta_s| = e^{-\pi\mu}$, has a saddle in the physical domain

$$ \tau_{\bullet}''(\beta) = \frac{-\mu}{\sqrt{k_{34}^2 - s^2}} \in C_I, \qquad (3.20) $$

while the saddle of the positive-frequency integral $F_R^{(3),\alpha}$ lies unphysically on the second Riemann sheet

$$ \tau_{\bullet}''(\alpha) = \frac{+\mu}{\sqrt{k_{34}^2 - s^2}} \in C_{II}. \qquad (3.21) $$

Consequently, we expect the contribution from both terms to be comparable, i.e. both are of order $O(e^{-\#\mu})$. Computing $F_R^{(3),\beta}$ is identical to that of $F_L^{(3),\alpha}$ and we directly write the result

$$ F_R^{(3),\beta} \approx \frac{\sqrt{\pi\mu} e^{-\pi\mu}}{(k_{34}^2 - s^2)^{3/4}} e^{-i\pi/4} \exp \left[ -i\mu \operatorname{arccosh} \frac{k_{34}}{s} \right] S(\tau_{\bullet}''(\beta)), \qquad (3.22) $$

where we have used the relation $\operatorname{Re} F_s(\tau_i) = \pi\mu$. As for the positive-frequency integral $F_R^{(3),\alpha}$, we need to continuously deform the integration contour so that it crosses the saddle point $\tau_\bullet''(\alpha)$ on the second Riemann sheet $C_{II}$. According to the saddle-point method, the complete integral is proportional to the integrand evaluated at the saddle point. Therefore, the size of the integral is estimated by

$$ |F_R^{(3),\alpha}| \sim \left| \exp \left( ik_{34} \tau_{\bullet}''(\alpha) + i \int_{-\infty}^{\tau_{\bullet}''(\alpha)} d\tau \, \omega_s(\tau) \right) \right| = \exp \left( - \operatorname{Im} \int_{-\infty}^{\tau_{\bullet}''(\alpha)} d\tau \, \omega_s(\tau) \right). \qquad (3.23) $$
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Analytic structure of the effective frequency $\omega_s$ in the complex $\tau''$ plane to compute $F_R^{(3)}$. For the positive-frequency integral $F_R^{(3),\alpha}$, the original integration contour is deformed to pass through the saddle point on the second Riemann sheet $\mathbb{C}_{II}$, which results in three different portions.</figcaption></figure>

The original integral contour is now deformed to pass through the saddle $\tau_{\bullet}''(\alpha) \in \mathbb{C}_{II}$. The new contour consists of three portions

$$ \int_{-\infty}^{\tau_{\bullet}''(\alpha)} = \int_{\mathbb{C}_1} + \int_{\mathbb{C}_2} + \int_{\mathbb{C}_3}, \qquad (3.24) $$

where the various contours are depicted in figure 9. Both $\mathbb{C}_1$ and $\mathbb{C}_3$ give purely real contributions, which directly drop out of our estimation. The only contribution to the imaginary part is from the semicircle $\mathbb{C}_2$, which can be nicely written as half the residue at the origin

$$ |F_R^{(3),\alpha}| \sim \exp \left( -\operatorname{Im} \int_{\mathbb{C}_2} d\tau \, \omega_s(\tau) \right) = e^{\pi \operatorname{Res}_{\tau=0} \omega_s(\tau)} = e^{-\pi \mu}. \qquad (3.25) $$

Interestingly, the price to pay for an “off-shell” vertex resonance is an exponential factor closely related to the location of the saddle point and the late-time behaviour of the effective frequency. We will see later that this is a general pattern. We can now finish the saddle-point approximation by performing the rest of the Gaussian integral. This yields

$$ F_R^{(3),\alpha} \approx \frac{-\sqrt{\pi\mu} e^{-\pi\mu}}{(k_{34}^2 - s^2)^{3/4}} e^{-i\pi/4} \exp \left[ -i\mu \operatorname{arccosh} \frac{k_{34}}{s} \right]. \qquad (3.26) $$

Combining (3.17), (3.22) and (3.26), and setting $S(\tau_{\bullet}''(\beta)) = 1$ in the soft limit, we obtain the full approximate solution for the cosmological collider signal
$$F_F + \text{c.c.} \approx \frac{2\pi\mu}{[(k_{12}^2 - s^2)(k_{34}^2 - s^2)]^{3/4}} \left\{ e^{-\pi\mu} \sin\mu \left[ \text{arccosh} \frac{k_{12}}{s} - \text{arccosh} \frac{k_{34}}{s} \right] + e^{-\pi\mu} \cos\mu \left[ \text{arccosh} \frac{k_{12}}{s} + \text{arccosh} \frac{k_{34}}{s} \right] \right\}, \quad (3.27)$$

for $k_{12} > k_{34}$. The opposite regime $k_{12} < k_{34}$ is obtained by swapping $k_{12} \leftrightarrow k_{34}$. This method gives very accurate results already for $\mu \gtrsim 2$, as we will see in section 3.3. The small phase shift for the non-local signal — oscillating in the $s^2/k_{12}k_{34}$ direction — can be fixed by considering higher-order corrections in $\mu^{-1}$.

**Discussion.** Let us make a few comments regarding the found approximate solution.

* The first line in (3.27) comes from the positive-frequency integral $F_R^{(3),\alpha}$ and asymptotes to $\sin[\mu \log(k_{12}/k_{34})]$ in the internal collapsed limit $s \to 0$. This term is known as the local cosmological collider signal as it is analytic in $s$. We recall that it captures the propagation of massive particles from one vertex to the other one [80]. More specifically, we see that the local signal is given by the time elapsed between the saddle point $\tau_{\bullet}''(\alpha)$ that creates the particle and the saddle $\tau_{\bullet}'(\alpha)$ that destroys the particle.

* The second line in (3.27) comes from the negative-frequency integral $F_R^{(3),\beta}$ and asymptotes to $\cos[\mu \log(4k_{12}k_{34}/s^2)]$ in the internal collapsed limit $s \to 0$. These oscillations are known as the non-local cosmological collider signal due to their non-analyticity in $s$. In contrast to the local signal, this one represents the propagation effect of massive particles from vacuum pair production to their subsequent decay at $\tau_{\bullet}'(\alpha)$ and $\tau_{\bullet}''(\beta)$ [80].

* In this minimal setup, both signals are suppressed by half the Boltzmann factor $e^{-\pi\mu} = e^{-m/(2T_{\text{ds}})}$, where $T_{\text{ds}} \equiv H/2\pi$ is the de Sitter temperature, but due to different reasons. For the local signal, it is due to the semi-circle we take to reach the saddle in the unphysical domain, whereas for the non-local signal, it is essentially the Bogoliubov coefficient. We will see in more general scenarios that the size of these signals can be drastically different.

* The seed integral we started with does not correspond to any physical coupling. However, all physical interaction vertices respecting the shift symmetry of external legs come with higher conformal time weights in the integrand, and their effects can be derived by differentiating the found result with respect to external kinematics. For instance, acting with $k_{12}\partial_{k_{12}}$ on (3.27) raises the conformal time power from zero to one. In the large-mass limit, the leading effect of such derivatives is to bring out powers of $\mu$ from the waveform while shifting its phase,

$$(k_{12}\partial_{k_{12}})^{p_1} (k_{34}\partial_{k_{34}})^{p_2} (F_F + \text{c.c.}) \approx \left( \frac{\mu k_{12}}{\sqrt{k_{12}^2 - s^2}} \right)^{p_1} \left( \frac{\mu k_{34}}{\sqrt{k_{34}^2 - s^2}} \right)^{p_2} \times (F_F + \text{c.c.})_{\text{shifted phase}}, \quad (3.28)$$
with

$$ \text{shifted phase} = \begin{cases} \frac{\pi(p_1 - p_2)}{2}, & \text{local signal} \\ \frac{\pi(p_1 + p_2)}{2}, & \text{non-local signal} \end{cases}. \qquad (3.29) $$

This feature does not alter the qualitative properties of the cosmological collider signals discussed above. In particular, the waveforms are unchanged modulo discrete phase shifts.

* Readers looking for mathematical rigour may have noticed that our “steepest” contour is, strictly speaking, invalid for the integral $F_R^{(3)}$. This is because the contour is restricted by the branching point $\tau_c^*$, and is in fact not a path for which the saddle-point approximation is justified. In other words, there exist points along this contour where the size of the integrand is actually larger than the height of the saddle at $\tau_\bullet''(\alpha)$. To remedy this seemingly fatal mathematical flaw, we note that the defining integrand given in terms of the exact mode function is analytic in the upper-half complex $\tau''$-plane without branching points. Moreover, the Bunch-Davies $i\epsilon$ prescription ensures that the integrand vanishes for the upper arc at infinity. Therefore, we can Wick rotate the contour from the very beginning, and then insert the WKB approximation. This way, the saddle can be picked up by the Wick rotated contour directly, eventually giving the same result as in (3.26).

### 3.2.2 Turning on a non-trivial sound speed

Let us now introduce a non-trivial sound speed. Recall that we have rescaled spacetime coordinates such that the speed of sound of the massive field remains unity while the external fields acquire a sound speed $c_s$, which is the ratio of both speeds. For $c_s > 1$, i.e. the exchanged massive field propagating slower than the external one, the analysis from the previous section applies and one only needs to replace $k_{12,34} \to c_s k_{12,34}$. The converse case $c_s < 1$, however, needs to be treated with care. Indeed, the bulk cutting rules require that both factorised integrals converge. Yet for $k_{12} > k_{34}$, convergence is achieved only if $k_{12}/s < c_s^{-1}$ as seen in section 2.2. In the opposite kinematic regime $k_{34} > k_{12}$, one needs to flip the other theta function, therefore requiring $k_{34}/s < c_s^{-1}$ for convergence. In combination, the convergence of the factorised time integrals for a reduced sound speed carves out a diamond-shape region in kinematic space to which our approximation scheme no longer applies, see the region in red in figure 10. Moreover, two new “edge” regions in green defined by $k_{34}/s < c_s^{-1} < k_{12}/s$ and $k_{12}/s < c_s^{-1} < k_{34}/s$ emerge where our approximate procedure still applies. From now on, we will assume $c_s < 1$ and focus on the edge region $k_{34}/s < c_s^{-1} < k_{12}/s$. The integral $F_L^{(3)}$ can again be separated into two pieces

$$ \begin{aligned} F_L^{(3)} & \approx \int_{-\infty+}^{0} \frac{d\tau'}{\sqrt{2\omega_s(\tau')}} \left[ e^{i(c_s k_{12}\tau' - \int_{\tau_i}^{\tau'} d\tau \omega_s(\tau))} + \beta_s(\tau')e^{i(c_s k_{12}\tau' + \int_{\tau_i}^{\tau'} d\tau \omega_s(\tau))} \right] \\ & \equiv F_L^{(3),\alpha} + F_L^{(3),\beta}, \end{aligned} \qquad (3.30) $$

Since the massive field remains untouched by the external reduced sound speed, the negative-frequency integral $F_L^{(3),\beta}$ remains subdominant and can be neglected. The saddle point for
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Kinematic space of the four-point function in the presence of a reduced sound speed for the external legs $c_s < 1$. The local signal direction is horizontal and the non-local signal direction is vertical.</figcaption></figure>

the positive-frequency integral $F_L^{(3),\alpha}$ is found to be

$$ \tau_{\bullet}^{\prime(\alpha)} = \frac{-\mu}{\sqrt{c_s^2 k_{12}^2 - s^2}} \in \mathbb{C}_I. \qquad (3.31) $$

Because $k_{12}/s > c_s^{-1}$, the saddle $\tau_{\bullet}^{\prime(\alpha)}$ lies in the physical region. Therefore, we can directly apply the result (3.17) after replacing $k_{12} \to c_s k_{12}$. Similarly, when $k_{34} > c_s^{-1}$, the second integral $F_R^{(3)}$ is also well approximated by (3.22) and (3.26) with $k_{34} \to c_s k_{34}$. However, in the edge region $k_{34}/s < c_s^{-1}$, the saddle points for the integral

$$ F_R^{(3)} \approx \int_{-\infty^+}^0 \frac{d\tau''}{\sqrt{2\omega_s(\tau'')}} \left[ e^{i(c_s k_{34} \tau'' + \int_{\tau_i}^{\tau''} d\tau \omega_s(\tau))} + \beta_s^*(\tau'') e^{i(c_s k_{34} \tau'' - \int_{\tau_i}^{\tau''} d\tau \omega_s(\tau))} \right] \quad (3.32) \\ \equiv F_R^{(3),\alpha} + F_R^{(3),\beta}, $$

lie along the imaginary axis

$$ \tau_{\bullet}''(\alpha) = \frac{+i\mu}{\sqrt{s^2 - c_s^2 k_{34}^2}} \in \mathbb{C}_{II}, \quad \tau_{\bullet}''(\beta) = \frac{+i\mu}{\sqrt{s^2 - c_s^2 k_{34}^2}} \in \mathbb{C}_I. \qquad (3.33) $$

Now, $\tau_{\bullet}''(\beta)$ is moved away from the physical regime. Using the rule of thumb, we expect $F_R^{(3),\beta}$ to be suppressed by $O(e^{-(\pi+\#)\mu})$, making it negligible compared to $F_R^{(3),\alpha} \sim O(e^{-\#\mu})$. We therefore drop the negative-frequency integral and focus on $F_R^{(3),\alpha}$. Following the saddle point approximation method, we deform the integration contour to pass through $\tau_{\bullet}''(\alpha)$. The
<figure><img src="image_15.png" /><figcaption>Figure 11. Analytic structure of the effective frequency $\omega_s$ in the complex $\tau''$ plane to compute $F_R^{(3)}$ in the case of a reduced sound speed $c_s < 1$ with $k_{34}/s < c_s^{-1}$.</figcaption></figure>

size of the integral is given by

$$ egin{aligned} |F_R^{(3),\alpha}| & \sim \left| \exp \left( i c_s k_{34} \tau_{\bullet}''(\alpha) + i \int_{-\infty}^{\tau_{\bullet}''(\alpha)} d\tau \, \omega_s(\tau) \right) \right| \\ & = \exp \left( -c_s k_{34} \operatorname{Im} \tau_{\bullet}''(\alpha) \right) \exp \left( -\operatorname{Im} \int_{-\infty}^{\tau_{\bullet}''(\alpha)} d\tau \, \omega_s(\tau) \right). \end{aligned} \quad (3.34) $$

To reach the imaginary saddle point on the second Riemann sheet, we need to rotate around the pole at the origin only a quarter of circle, plus an imaginary contour from the branching point to the saddle point. The new contour consists now of four portions

$$ \int_{-\infty}^{\tau_{\bullet}''(\alpha)} = \int_{C_1} + \int_{C_2} + \int_{C_3} + \int_{C_4}, \quad (3.35) $$

where the various contours are depicted in figure 11. Picking up a quarter of residue at the origin, we obtain

$$ egin{aligned} |F_R^{(3)}| & \sim \exp \left( -\operatorname{Im} \int_{C_2} d\tau \, \omega_s(\tau) \right) \exp \left( -\operatorname{Im} \int_{C_4} d\tau \, \omega_s(\tau) \right) \exp \left( -c_s k_{34} \operatorname{Im} \tau_{\bullet}''(\alpha) \right) \\ & = e^{\frac{\pi}{2} \operatorname{Res}_{\tau=0} \omega_s(\tau)} \exp \left( \int_{\mu/s}^{\mu/\sqrt{s^2 - c_s^2 k_{34}^2}} \sqrt{s^2 - \frac{\mu^2}{x^2}} dx - \frac{c_s k_{34} \mu}{\sqrt{s^2 - c_s^2 k_{34}^2}} \right) \quad (3.36) \\ & = e^{-\pi \mu/2} \exp \left( -\mu \arcsin \frac{c_s k_{34}}{s} \right). \end{aligned} $$

Applying the full saddle point approximation yields

$$ F_R^{(3)} = \frac{-i \sqrt{\pi \mu} e^{-\pi \mu/2}}{(s^2 - c_s^2 k_{34}^2)^{3/4}} \exp \left[ -\mu \arcsin \frac{c_s k_{34}}{s} \right]. \quad (3.37) $$
Combining all relevant expressions, we obtain the full approximate solution for the cosmological collider signal

$$F_F + \text{c.c.} \approx \frac{2\pi\mu}{[(c_s^2 k_{12}^2 - s^2)(s^2 - c_s^2 k_{34}^2)]^{3/4}} \times e^{-\pi\mu/2} \exp\left[-\mu \arcsin\left(\frac{c_s k_{34}}{s}\right)\right] \cos\left[\mu \operatorname{arccosh}\left(\frac{c_s k_{12}}{s} + \frac{3\pi}{4}\right)\right]. \quad (3.38)$$

Again, the signal in the edge region $k_{12}/s < c_s^{-1} < k_{34}/s$ is obtained by swapping the kinematic variables $k_{12} \leftrightarrow k_{34}$.

**Discussion.** A few comments are in order.

* In contrast to the conventional cosmological collider signals, the one with a reduced sound speed is suppressed by only a *quarter* of the usual Boltzmann factor $e^{-\pi\mu/2} = e^{-m/(4T_{ds})}$ due to the *quarter* circle we take to reach the saddle. Such amplified oscillations have long been known in the literature but without much explanation nor interpretation [11, 43].

* The behaviour of the cosmological collider signal is intimately related to the location of the saddle $\tau_{\bullet}^{(\alpha)}$. In the edge region of kinematic space, $\tau_{\bullet}^{(\alpha)}$ is real, and it shifts in the real direction when we vary the ratio $k_{12}/s$, giving rise to oscillations. On the other hand, $\tau_{\bullet}^{(\alpha)}$ is purely imaginary. Hence when we vary the ratio $k_{34}/s$, the saddle point shifts in the imaginary direction, giving rise to smooth exponential attenuation instead of oscillations. In addition, the arcsin function nicely produces a smooth interpolation between $e^{-\pi\mu/2}$ (reduced sound speed) and $e^{-\pi\mu}$ (conventional collider), as $c_s k_{34}/s$ increases from zero to one.

* Interestingly, we notice that (3.38) can be derived from (3.27) via direct analytic continuation. To see this, note that the arccosh ($z$) function has a branching point at $z \equiv c_s k_{34}/s = 1$. To consistently continue beyond this singularity, we recall that $k_{34} = |k_{34}| - i\epsilon$ is defined with a negative imaginary part to ensure the convergence of $F_R^{(3)}$. Therefore, the correct transformation of the arccosh ($z$) function is

$$\operatorname{arccosh}(z) \rightarrow -i \operatorname{arccos}(z), \quad 0 < \operatorname{Re} z < 1, \quad \operatorname{Im} z < 0. \quad (3.39)$$

The power law prefactor transforms as

$$(c_s^2 k_{34}^2 - s^2)^{3/4} \rightarrow e^{-3i\pi/4} (s^2 - c_s^2 k_{34}^2)^{3/4}. \quad (3.40)$$

Expansion of the trigonometric functions yields terms of higher order in $O(e^{-\#\mu})$ which can be dropped. This allows us to directly land on the correct result (3.38) without repeating the entire saddle point analysis. It is intriguing why this should be the case since approximate solutions are not expected to enjoy precise analytic continuations in general. We leave this pleasant yet seemingly coincidental phenomenon for future work.

* Similar to the case with unit sound speed, realistic correlators with more derivatives in the interaction vertices can be obtained by differentiating the seed integral with respect
to the external kinematics. In the large-$\mu$ asymptotic limit, we can approximate

$$ (k_{12} \partial_{k_{12}})^{p_1} (k_{34} \partial_{k_{34}})^{p_2} (F_F + \text{c.c.}) \approx \left( \frac{c_s \mu k_{12}}{\sqrt{c_s^2 k_{12}^2 - s^2}} \right)^{p_1} \left( \frac{-c_s \mu k_{34}}{\sqrt{s^2 - c_s^2 k_{34}^2}} \right)^{p_2} \times (F_F + \text{c.c.})_{\text{shifted phase}}. \quad (3.41) $$

Again, the waveform changes at most by a discrete phase shift,

$$ \text{shifted phase} = \frac{\pi p_1}{2}. \quad (3.42) $$

Consequently, the phase of the low-speed collider signal is fixed at $\pm\pi/4$, $\pm3\pi/4$ at leading order in $\mu^{-1}$, which is a unique prediction of the approximated signal.

### 3.2.3 Further turning on the chemical potential

We now move to the most general case where we further turn on a non-zero chemical potential $\kappa$. Since the WKB approximation (3.4) is valid only when $|\beta_s(\tau)| \ll 1$, we will only consider the parameter regime with weak particle production, i.e. $\kappa < \mu$. This avoids large backreaction from $\sigma$ loops [71]. Without loss of generality, we set $k_{12} > k_{34}$ and concentrate on the enhanced helicity $\lambda = -1$. The result for the $\lambda = +1$ mode can be found via a simple flip of the chemical potential $\kappa \to -\kappa$. As discussed in the previous subsection, we also require $k_{12}/s > c_s^{-1}$ to ensure the validity of the bulk cutting rule. The effect of the chemical potential is to add a linear term in the massive effective frequency

$$ \omega_s^2(\tau) = s^2 + \frac{2s\kappa}{\tau} + \frac{\mu^2}{\tau^2}. \quad (3.43) $$

This additional linear term does not alter the leading late-time pole at $\tau = 0$, but it shifts the turning points away from the imaginary axis,

$$ \tau_c = \frac{1}{s} \left( -\kappa - i \sqrt{\mu^2 - \kappa} \right). \quad (3.44) $$

The branch cut at $\omega_s^2(\tau) < 0$ joins $\tau_c, \tau_c^*$ and 0, separating the complex plane into different Riemann sheets. The effective frequency is therefore defined as double-valued function

$$ \begin{aligned} \omega_s(\tau) &= +\sqrt{s^2 + \frac{2s\kappa}{\tau} + \frac{\mu^2}{\tau^2}}, & \tau \in \mathbb{C}_I, \\ \omega_s(\tau) &= -\sqrt{s^2 + \frac{2s\kappa}{\tau} + \frac{\mu^2}{\tau^2}}, & \tau \in \mathbb{C}_II. \end{aligned} \quad (3.45) $$

The singulant function admits a closed-form expression

$$ F_s(\tau) = -2i \left[ \tau\omega_s(\tau) + \mu \operatorname{arctanh} \left( \frac{s\tau\kappa + \mu^2}{-\mu\tau\omega_s(\tau)} \right) - \kappa \operatorname{arctanh} \left( \frac{s\tau + \kappa}{-\tau\omega_s(\tau)} \right) \right] + \pi(\mu - \kappa). \quad (3.46) $$
The chemical potential assists spontaneous particle production by alleviating the effective mass in the Boltzmann factor $|β_s|^2 = e^{-2\text{Re}F_s(τ_i)} = e^{-2π(μ-κ)}$. Following the same procedure as for the previous cases, we start with the left time integral

$$F_L^{(3)} \approx \int_{-\infty^+}^0 \frac{d\tau'}{\sqrt{2\omega_s(\tau')}} \left[ e^{i(k_{12}\tau' - \int_{\tau_i}^{\tau'} d\tau \omega_s(\tau))} + \beta_s(\tau')e^{i(k_{12}\tau' + \int_{\tau_i}^{\tau'} d\tau \omega_s(\tau))} \right] \quad (3.47) \\ \equiv F_L^{(3),\alpha} + F_L^{(3),\beta}.$$

Solving for the saddle points, we find

$$\tau_{\bullet}^{\prime(\alpha)} = \frac{s(\kappa - \mu\Delta_L)}{c_s^2 k_{12}^2 - s^2} \in \mathbb{C}_I, \quad \tau_{\bullet}^{\prime(\beta)} = \frac{s(\kappa + \mu\Delta_L)}{c_s^2 k_{12}^2 - s^2} \in \mathbb{C}_II, \quad (3.48)$$

with $\Delta_L = \sqrt{c_s^2 k_{12}^2/s^2 - (1 - \kappa^2/\mu^2)}$. Again, the integral $F_L^{(3),\alpha}$ has a physical saddle point, and thereby dominates over $F_L^{(3),\beta}$. Finishing the Gaussian integral around $\tau_{\bullet}^{\prime(\alpha)}$, we obtain

$$F_L^{(3)} \approx \frac{\sqrt{\pi}\mu^{3/2}}{s^{3/2}(\kappa + \mu\Delta_L)\sqrt{\Delta_L}} e^{-i\pi/4 - i\theta} \exp \left[ -i\mu \left( \text{arccosh} \frac{c_s k_{12}/s}{\sqrt{1 - \kappa^2/\mu^2}} \right) -i\kappa \left( \text{arccoth} \frac{c_s k_{12}}{s} - \text{arctanh} \frac{c_s k_{12}\kappa}{s\mu\Delta_L} \right) \right], \quad (3.49)$$

where $\theta$ is a kinematic-independent phase, given by

$$\theta(\mu, \kappa) \equiv \frac{\mu + \kappa}{2} \log(\mu + \kappa) - \frac{\mu - \kappa}{2} \log(\mu - \kappa) - \kappa. \quad (3.50)$$

As usual, the integral $F_R^{(3)}$ is more complicated. We define

$$F_R^{(3)} \approx \int_{-\infty^+}^0 \frac{d\tau''}{\sqrt{2\omega_s(\tau'')}} \left[ e^{i(c_s k_{34}\tau'' + \int_{\tau_i}^{\tau''} d\tau \omega_s(\tau))} + \beta_s^*(\tau'')e^{i(c_s k_{34}\tau'' - \int_{\tau_i}^{\tau''} d\tau \omega_s(\tau))} \right] \quad (3.51) \\ \equiv F_R^{(3),\alpha} + F_R^{(3),\beta}.$$

The location of the saddle points are

$$\tau_{\bullet}^{\prime\prime(\alpha)} = \frac{s(\kappa + \mu\Delta_R)}{c_s^2 k_{34}^2 - s^2} \in \mathbb{C}_II, \quad \tau_{\bullet}^{\prime\prime(\beta)} = \frac{s(\kappa - \mu\Delta_R)}{c_s^2 k_{34}^2 - s^2} \in \mathbb{C}_I, \quad (3.52)$$

with $\Delta_R = \sqrt{c_s^2 k_{34}^2/s^2 - (1 - \kappa^2/\mu^2)}$. For $k_{34}/s > c_s^{-1}\sqrt{1 - \kappa^2/\mu^2}$, we have a real $\Delta_R$ and the saddle points lie on the real axis. For $k_{34}/s < c_s^{-1}\sqrt{1 - \kappa^2/\mu^2}$, $\Delta_R$ is purely imaginary and the saddle points can move around in the complex plane. The point $k_{34}/s = c_s^{-1}$ also distinguishes the position of $\tau_{\bullet}^{\prime\prime(\alpha)}$ relative to infinity. As the qualitative behaviour of cosmological collider signals dramatically depends on the location of the saddle points, we separate the discussion according to kinematic regions.

**Interior region $k_{34}/s > c_s^{-1}$**. This region is depicted in blue in figure 10 (recall that $k_{12} > k_{34}$). In this region, $\tau_{\bullet}^{\prime\prime(\alpha)}$ is along the positive real axis on the second Riemann sheet
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Analytic structure of the effective frequency $\omega_s$ in the complex $\tau''$ plane to compute $F_R^{(3)}$ in the case of a non-unity sound speed and non-zero chemical potential in the interior region.</figcaption></figure>

$C_{II}$, whereas $\tau_{\bullet}''^{(\beta)}$ is along the negative real axis on $C_{I}$. The integral $F_R^{(3),\beta}$ picks up the saddle in the physical region, which leads to

$$F_R^{(3),\beta} \approx \frac{\sqrt{\pi}\mu^{3/2}e^{-\pi(\mu-\kappa)}}{s^{3/2}(\kappa + \mu\Delta_R)\sqrt{\Delta_R}} e^{i\pi/4+i\theta} \exp \left[ -i\mu \left( \operatorname{arccosh} \frac{c_s k_{34}/s}{\sqrt{1-\kappa^2/\mu^2}} \right) -i\kappa \left( \operatorname{arccoth} \frac{c_s k_{34}}{s} - \operatorname{arctanh} \frac{c_s k_{34}\kappa}{s\mu\Delta_R} \right) \right]. \quad (3.53)$$

In order to apply the saddle-point approximation for the integral $F_R^{(3),\alpha}$, we deform the initial integration contour such that it dives under the branch cut, passing through $\tau_{\bullet}''^{(\alpha)}$ in the steepest direction. As represented in figure 12, we again need to rotate around the origin, therefore picking up the same suppression power

$$|F_R^{(3),\alpha}| \sim \exp \left( -\operatorname{Im} \int_{C_2} d\tau \, \omega_s(\tau) \right) = e^{\frac{\pi}{\operatorname{Res}_{\tau=0}} \omega_s(\tau)} = e^{-\pi\mu}. \quad (3.54)$$

Finishing the Gaussian integral, we obtain

$$F_R^{(3),\alpha} \approx \frac{-\sqrt{\pi}\mu^{3/2}e^{-\pi\mu}}{s^{3/2}(\kappa - \mu\Delta_R)\sqrt{\Delta_R}} e^{i\pi/4+i\theta} \exp \left[ i\mu \left( \operatorname{arccosh} \frac{c_s k_{34}/s}{\sqrt{1-\kappa^2/\mu^2}} \right) -i\kappa \left( \operatorname{arccoth} \frac{c_s k_{34}}{s} + \operatorname{arctanh} \frac{c_s k_{34}\kappa}{s\mu\Delta_R} \right) \right]. \quad (3.55)$$

Combining all found results, we obtain the leading approximation of the local and non-local signals in the interior region $k_{12}/s > k_{34}/s > c_s^{-1}$,
$$F_F + \text{c.c.} \approx \frac{2\pi\mu^3}{s^3(\kappa + \mu\Delta_L)\sqrt{\Delta_L\Delta_R}} \times \left\{ \frac{e^{-\pi\mu}}{\mu\Delta_R - \kappa} \sin\left[ \mu \left( \operatorname{arccosh} \frac{c_s k_{12}/s}{\sqrt{1-\kappa^2/\mu^2}} - \operatorname{arccosh} \frac{c_s k_{34}/s}{\sqrt{1-\kappa^2/\mu^2}} \right) + \Xi \right] + \frac{e^{-\pi(\mu-\kappa)}}{\mu\Delta_R + \kappa} \cos\left[ \mu \left( \operatorname{arccosh} \frac{c_s k_{12}/s}{\sqrt{1-\kappa^2/\mu^2}} + \operatorname{arccosh} \frac{c_s k_{34}/s}{\sqrt{1-\kappa^2/\mu^2}} \right) + \Upsilon \right] \right\}, \quad (3.56)$$

where the dephasing functions are given by

$$\Xi \equiv \kappa \left( \operatorname{arccoth} \frac{c_s k_{12}}{s} + \operatorname{arccoth} \frac{c_s k_{34}}{s} - \operatorname{arctanh} \frac{c_s k_{12}\kappa}{s\mu\Delta_L} + \operatorname{arctanh} \frac{c_s k_{34}\kappa}{s\mu\Delta_R} \right), \quad (3.57) \\ \Upsilon \equiv \Xi - 2\kappa \operatorname{arctanh} \frac{c_s k_{34}\kappa}{s\mu\Delta_R}.$$

Let us now comment on the result.

* In the interior region, the cosmological collider signals here behave qualitatively the same way as those in the case $c_s = 1$ and $\kappa = 0$. The first line in (3.56) represents the local signal. Its size is dictated only by the mass $O(e^{-\pi\mu})$ coming from the semicircle around the origin. In the soft limit, the dephasing function approaches zero $\Xi \to 0$ and the waveform approaches $\sin\mu\log(k_{12}/k_{34})$. In other words, the phase of the local cosmological collider signal is universally zero, independent of any parameters of the theory.

* The second line in (3.56) gives the non-local signal, whose size is amplified by the chemical potential $O(e^{-\pi(\mu-\kappa)})$. The corresponding dephasing function asymptotes to a constant that does depend on the parameters of the theory $\Upsilon \to -2\kappa \operatorname{arctanh} \kappa/\mu$. In general, the presence of the dephasing function implies that the waveform can significantly deviate from a simple logarithmically oscillating pattern.

* As in the previous cases, results from realistic interaction vertices are straightforwardly obtained by taking derivatives with respect to external kinematics. The waveforms stay invariant modulo discrete phase shifts.

**Edge region $k_{34}/s < c_s^{-1}\sqrt{1-\kappa^2/\mu^2}$**. This region is depicted in green in figure 10. Here, both saddle points are complex. Using the rule of thumb, the saddle at $\tau_{\bullet}''(\beta)$ can then be dropped due to its additional suppression. We then focus on the saddle located at

$$\tau_{\bullet}''(\alpha) = \frac{s(-\kappa + i\mu|\Delta_R|)}{s^2 - c_s^2 k_{34}^2} \in \mathbb{C}_{\text{II}}. \quad (3.58)$$
To reach this saddle, we again need to dive under the branch cut to enter the second Riemann sheet $\mathcal{C}_{II}$. The size of the integral is estimated from the exponent

$$ egin{aligned} |F_R^{(3),\alpha}| & \sim \left| \exp \left( i c_s k_{34} \tau_{\bullet}''^{(\alpha)} + i \int_{-\infty}^{\tau_{\bullet}''^{(\alpha)}} d\tau \, \omega_s(\tau) \right) \right| \\ & = \exp \left( -c_s k_{34} \operatorname{Im} \tau_{\bullet}''^{(\alpha)} \right) \exp \left( -\operatorname{Im} \int_{-\infty}^{\tau_{\bullet}''^{(\alpha)}} d\tau \, \omega_s(\tau) \right). \end{aligned} \quad (3.59) $$

As depicted in figure 13, the path to reach the saddle consists of four segments

$$ \int_{-\infty}^{\tau_{\bullet}''^{(\alpha)}} = \int_{C_1} + \int_{C_2} + \int_{C_3} + \int_{C_4}. \quad (3.60) $$

The contour $C_1$ is purely real and drops out of the imaginary part. The contour $C_4$ part is negligible if $k_{34}/s \ll c_s^{-1}$ because both the integrand and the length of the integration path tend to zero in this regime. We are left with an integral along a quarter circle around the origin, which gives

$$ \operatorname{Im} \int_{C_2} d\tau \, \omega_s(\tau) = - \operatorname{Im} \left( \frac{i\pi}{2} \operatorname{Res}_{\tau=0} \omega_s(\tau) \right) = \frac{\pi\mu}{2}, \quad (3.61) $$

and a contour around the branch cut on $\mathcal{C}_{II}$, whose imaginary part can be computed as

$$ \begin{aligned} \operatorname{Im} \int_{C_3} d\tau \, \omega_s(\tau) &= \frac{1}{4} \operatorname{Im} \int_{C_{\text{cut}} \subset \mathcal{C}_{II}} d\tau \, \omega_s(\tau) = \frac{1}{4} \operatorname{Im} \int_{C_{\infty} \subset \mathcal{C}_{II}} d\tau \, \omega_s(\tau) \\ &= \operatorname{Im} \left( \frac{-\pi i}{2} \operatorname{Res}_{\tau=\infty} \omega_s(\tau) \right) = -\frac{\pi \kappa}{2}. \end{aligned} \quad (3.62) $$

Here in the first step we have used the fact that $\omega_s^*(\tau) = \omega_s(\tau)$ is a real function and that along the two sides of the branch cut, $\omega_s(\tau)$ differs by a minus sign. In the second step, we have deformed the contour around the branch cut $C_{\text{cut}}$ to a large circle at infinity $\mathcal{C}_{\infty}$ on Sheet II, which is then evaluated as a residue at infinity. We thus obtain a neat size formula $|F_R^{(3),\alpha}| \sim e^{-\pi(\mu-\kappa)/2}$. We note in passing that it is remarkable to behold the intimate connection between the cosmological collider signal strength, which reflects the amount of kinematic non-analyticities in boundary correlators, and the residue of the dispersion relation at early and late times, which reflects the amount of temporal non-analyticities of the bulk Lagrangian. Completing the Gaussian integral, we then have

$$ F_R^{(\alpha)} \approx \frac{\sqrt{\pi} \mu^{3/2} e^{-\pi(\mu-\kappa)/2}}{s^{3/2}(\kappa + i\mu|\Delta_R|)\sqrt{|\Delta_R|}} e^{i\theta} \exp \left[ -\mu \arcsin \frac{c_s k_{34}/s}{\sqrt{1-\kappa^2/\mu^2}} \right. \\ \left. + \kappa \arctan \frac{c_s k_{34}\kappa}{s\mu|\Delta_R|} - i\kappa \operatorname{arctanh} \frac{c_s k_{34}}{s} \right]. \quad (3.63) $$

Combining all found results, we obtain the leading approximation of the signals in the edge region
<figure><img src="image_17.png" /><figcaption>Figure 13. Analytic structure of the effective frequency $\omega_s$ in the complex $\tau''$ plane to compute $F_R^{(3)}$ in the case of a non-unity sound speed and non-zero chemical potential in the edge region.</figcaption></figure>

$$ F_F + \text{c.c.} \approx \frac{2\pi\mu^3 e^{-\pi(\mu-\kappa)/2}}{s^3(\kappa + \mu\Delta_L)[(\kappa^2 + \mu^2|\Delta_R|^2)\Delta_L|\Delta_R|]^{1/2}} \\ \times \exp\left[-\mu \arcsin\frac{c_s k_{34}/s}{\sqrt{1-\kappa^2/\mu^2}} + \kappa \arctan\frac{c_s k_{34}\kappa}{s\mu|\Delta_R|}\right] \\ \times \cos\left[\mu \operatorname{arccosh}\frac{c_s k_{12}/s}{\sqrt{1-\kappa^2/\mu^2}} + \kappa \operatorname{arctanh}\frac{c_s k_{34}}{s} + \frac{3\pi}{4} + \Omega\right], \quad (3.64) $$

where the dephasing function reads

$$ \Omega \equiv \kappa \left( \operatorname{arccoth} \frac{c_s k_{12}}{s} - \operatorname{arctanh} \frac{c_s k_{12} \kappa}{s \mu \Delta_L} \right) - \operatorname{arctan} \frac{\kappa}{\mu |\Delta_R|}. \quad (3.65) $$

* In the edge region of the kinematics space, the cosmological collider signal has important chemical-potential-dependent corrections. Apart from the overall enhancement, the corrected attenuation factor renders the signal amplitude to smoothly interpolate between $e^{-\pi(\mu-\kappa)/2}$ and $e^{-\pi(\mu-\kappa)}$ when $c_s k_{34}/s$ varies from 0 to $\sqrt{1-\kappa^2/\mu^2}$. This corresponds to the shift of the saddle $\tau_{\bullet}''(\alpha)$ in the imaginary direction.

* The waveform is also corrected by the dephasing function $\Omega$, which tends to a constant in the hierarchical soft limit $c_s k_{12}/s \to \infty$, $c_s k_{34}/s \to 0$,

$$ \lim_{\substack{c_s k_{12}/s \to \infty \\ c_s k_{34}/s \to 0}} \Omega = -\kappa \operatorname{arctanh} \frac{\kappa}{\mu} - \operatorname{arctan} \frac{\kappa}{\sqrt{1 - \kappa^2/\mu^2}}. \quad (3.66) $$

Interestingly, in addition to the dephasing function, the waveform acquires a term proportional to the chemical potential (the second term in the second line of (3.64)).
With $k_{12}$ held fixed, there is an intermediate stage where $k_{34}/s \ll c_s^{-1}$ and the waveform is approximated by

$$F_F + \text{c.c.} \supset \exp\left(-c_s\sqrt{\mu^2 - \kappa^2}\frac{k_{34}}{s} + \dots\right) \cos\left(c_s\kappa\frac{k_{34}}{s} + \dots\right). \quad (3.67)$$

Consequently, a new type of cosmological collider signal arises. Drastically distinct from all other more standard signals above, this new signal oscillates linearly with respect to the momentum ratio $k_{34}/s$ with a frequency $c_s\kappa$ determined by the sound speed and the chemical potential. Meanwhile, the amplitude of the oscillations decays exponentially at a rate $c_s\sqrt{\mu^2 - \kappa^2}$ which is often comparable to oscillation frequency. Due to its large decay rate, we name this new signal the transient cosmological collider signal. Such a transient signal was found previously in the parity-odd sector of non-local EFT description of the same system [27], but without apparent physical interpretation. Here we see clearly that the origin of the transient signal is the curvilinear motion of the saddle point $\tau_{\bullet}''(\alpha)$ on the complex plane. When the ratio $c_s k_{34}/s$ increases from nearly zero, the saddle point $\tau_{\bullet}''(\alpha)$ moves away from $\tau_c^*$ and shifts simultaneously in the real and imaginary direction, leading to the damped oscillations in the transient signal. The oscillations stop exponentially decaying when $\tau_{\bullet}''(\alpha)$ hits the negative real axis on the second Riemann sheet $\mathbb{C}_II$, which corresponds to the boundary of the edge region where $k_{34}/s = c_s^{-1}\sqrt{1 - \kappa^2/\mu^2}$ and our approximation fails.

* Similar to the previous case, the found result (3.64) can be derived from direct analytic continuation from that in the interior region (3.56), upon utilising the following continuation rules,

$$\begin{align}
\operatorname{arccosh} \frac{c_s k_{34}/s}{\sqrt{1 - \kappa^2/\mu^2}} &\rightarrow -i \operatorname{arccos} \frac{c_s k_{34}/s}{\sqrt{1 - \kappa^2/\mu^2}}, \nonumber \\
\Delta_R &\rightarrow -i |\Delta_R|, \tag{3.68} \\
\operatorname{arccoth} \frac{c_s k_{34}}{s} &\rightarrow \operatorname{arctanh} \frac{c_s k_{34}}{s} + \frac{i\pi}{2}. \nonumber
\end{align}$$

**Joint region $c_s^{-1}\sqrt{1 - \kappa^2/\mu^2} < k_{34}/s < c_s^{-1}$**. Entering the join region from above, $k_{34} > s/c_s \rightarrow k_{34} < s/c_s$, the location of $\tau_{\bullet}''(\alpha)$ switches from $+\infty$ to $-\infty$, while the location of $\tau_{\bullet}''(\beta)$ changes continuously. In principle, we can run the same reasoning as before and perform the same saddle-point analysis. To avoid technical repetitions, instead of computing the approximate results from scratch, we shall take advantage of the analyticity property we discovered above. The continued results are then checked against the exact results. We observe that the dephasing functions contain singularities at $k_{34}/s = c_s^{-1}$ (where $\Delta_R = \tilde{\kappa}/\mu$) for different terms. In $\Upsilon$, the singularities of two trigonometrical functions cancel out, leaving a straightforward analytic continuation rule

$$\Upsilon: \left\{ \begin{array}{l} \operatorname{arccoth} \frac{c_s k_{34}}{s} \rightarrow \operatorname{arctanh} \frac{c_s k_{34}}{s} + \frac{i\pi}{2} \\ \\ \operatorname{arctanh} \frac{c_s k_{34}\kappa}{s\mu\Delta_R} \rightarrow \operatorname{arccoth} \frac{c_s k_{34}\kappa}{s\mu\Delta_R} + \frac{i\pi}{2} \end{array} \right. \quad (3.69)$$
In $\Xi$, the singularities add up, suggesting that the consistent continuation rule is

$$ \Xi: \left\{ \begin{array}{l} \text{arccoth} \frac{c_s k_{34}}{s} \rightarrow \text{arctanh} \frac{c_s k_{34}}{s} + \frac{i\pi}{2} \\ \\ \text{arctanh} \frac{c_s k_{34} \kappa}{s\mu\Delta_R} \rightarrow \text{arccoth} \frac{c_s k_{34} \kappa}{s\mu\Delta_R} - \frac{i\pi}{2} \end{array} \right. \qquad (3.70) $$

All other terms remain unchanged.

## 3.3 Refined cosmological collider templates

Unravelling the rich bulk physics hidden in the boundary correlators requires searching for cosmological collider signals as intricate kinematic dependences in the non-Gaussian statistics of cosmological perturbations. Due to the high dimensionality of the CMB/LSS datasets, reconstructing the correlators in the full kinematic domain using an unbiased estimator is extremely challenging. Instead, one is compelled to use biased estimators which are derived or motivated from specific theoretical models. In other words, this amounts to computing the inner product between a shape template and the data to check for the detection significance, which in turn requires an efficient sampling of the high-dimensional kinematics space as well as parameter space. As seen in section 2, exact solutions of the cosmological collider signals often involve sums of products of hypergeometric functions, leading to significant computational complexity. This motivates us to propose approximate shape templates that substantially reduce computational cost while maintaining accuracy across most of the relevant kinematic and parameter spaces.

The saddle point method has proved to be a powerful approximation scheme that allowed us to link kinematic dependence of the boundary correlators to the resonance history in the inflationary bulk. As a useful by-product, it has provided clean boost-breaking templates expressed in terms of elementary functions, yet more accurate than the ones previously used. Here, we propose new refined templates for cosmological collider signals that match exact calculations both in terms of amplitude and function dependence of the waveform. Importantly, as probing ultra-soft limits of cosmological correlators requires a large hierarchy of scales and is limited by cosmic variance, these new templates are valid in extended and intermediate regions of the parameters space (e.g. covering mildly reduced sound speed cases) and of the kinematic configurations (e.g. covering mildly soft limits).

**Trispectrum shape templates.** Following standard conventions, we define the dimensionless trispectrum shape function by

$$ \langle \zeta_{k_1} \zeta_{k_2} \zeta_{k_3} \zeta_{k_4} \rangle \equiv (2\pi)^9 \delta^3(\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3 + \mathbf{k}_4) \frac{(k_T/4)^3 P_{\zeta}^3}{(k_1 k_2 k_3 k_4)^3} \mathcal{T}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4), \qquad (3.71) $$

where $k_T \equiv \sum_i k_i$ and $P_\zeta \equiv \frac{k^3}{2\pi^2} P_\zeta \approx 2.1 \times 10^{-9}$ is the dimensionless power spectrum of the curvature perturbation $\zeta$. The function $\mathcal{T}$ is a dimensionless scalar function of momenta that is invariant under the (same) rescaling of all momenta $\mathbf{k}_i \rightarrow \alpha \mathbf{k}_i$. The dimensionless trispectrum $\mathcal{T}$ can be further decomposed into a polarisation factor $\Pi$ and a waveform function $\mathcal{F}$

$$ \mathcal{T}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4) = \left(\frac{s^2}{k_{12}k_{34}}\right)^{3/2} \sum_{\lambda} \Pi_{S,\lambda}(\hat{\mathbf{k}}_1, \hat{\mathbf{k}}_2, \hat{\mathbf{k}}_3, \hat{\mathbf{k}}_4, \hat{s}) \mathcal{F}^{\lambda}(k_{12}, k_{34}, s) + 2 \text{ perms}, \quad (3.72) $$
where $\hat{k}_i \equiv k_i/k_i$, $\hat{s} \equiv s/s$ are unit vectors, and $s \equiv k_1 + k_2$. The two permutations correspond to the t- and u-channels, with $t \equiv k_1 + k_3$ and $u \equiv k_1 + k_4$. In the following, we only consider the dominant helicity component of $\mathcal{F}^\lambda$. The polarisation factor $\Pi_{S,\lambda}$ depends on the (integer) spin of the exchanged field $S$, the helicity $-S \le \lambda \le S$ and the detailed form of the interaction vertices. In the soft limit, the polarisation factor is well approximated by [7]

$$ \Pi_{S,\lambda}(\hat{k}_1, \hat{k}_2, \hat{k}_3, \hat{k}_4, \hat{s}) \approx e^{i\lambda\psi} P_S^\lambda(\hat{k}_1 \cdot \hat{s}) P_S^{-\lambda}(\hat{k}_3 \cdot \hat{s}), \quad (3.73) $$

where $P_S^\lambda$ is the associated Legendre polynomial, and $\psi$ is the dihedral angle between the planes spanned by $k_1, k_2$ and $k_3, k_4$, such that $\cos \psi \equiv \hat{n}_{12} \cdot \hat{n}_{34}$ with $\hat{n}_{12} \equiv \hat{k}_1 \times \hat{k}_2$ (resp. $\hat{n}_{34} \equiv \hat{k}_3 \times \hat{k}_4$) being the normal unit vector to the plane generated by the vectors $k_1, k_2$ (resp. $k_3, k_4$). Explicitly, the first few terms are $\Pi_{0,0} = 1, \Pi_{1,0} = (\hat{k}_1 \cdot \hat{s})(\hat{k}_3 \cdot \hat{s})$, and $\Pi_{1,\pm 1} = -\frac{1}{2}e^{\pm i\psi}\sqrt{[1-(\hat{k}_1 \cdot \hat{s})^2][1-(\hat{k}_3 \cdot \hat{s})^2]}$. The s-channel dimensionless waveform function $\mathcal{F}(k_{12}, k_{34}, s)$ can be decomposed into a local part that exhibits an oscillatory behaviour in the logarithm of the momentum ratio $k_{12}/k_{34}$ (assuming $k_{12} > k_{34}$) and a non-local part oscillating in the logarithm of the momentum ratio $k_{12}k_{34}/s^2$:<fnref n="10" />

$$ \mathcal{F}(k_{12}, k_{34}, s) = \mathcal{F}^{\mathrm{L}}(k_{12}, k_{34}, s) + \mathcal{F}^{\mathrm{NL}}(k_{12}, k_{34}, s). \quad (3.74) $$

The waveforms $\mathcal{F}^\mathrm{L}$ and $\mathcal{F}^\mathrm{NL}$ are parametrised by three constants: (i) a dimensionless mass parameter $\mu$ of the exchanged massive field, (ii) a speed of sound<fnref n="11" /> $0 < c_s < \infty$, and (iii) a dimensionless chemical potential $|\kappa| < \mu$. In the s-channel, all cosmological collider signals lie in the soft limit $s \to 0$. Without loss of generality, we set $k_{12} > k_{34}$ and the opposite case is obtained by the trivial swap $k_{12} \leftrightarrow k_{34}$. With the presence of non-trivial sound speed $c_s \ne 1$ and a non-zero chemical potential $\kappa \ne 0$, which signals the breaking of de Sitter boosts, the kinematics space can be divided into an interior region $k_{34}/s \ge c_s^{-1}$ and an edge region $k_{34}/s \le c_s^{-1}$, see figure 10. The template then takes the form of a piecewise function of $k_{12}, k_{34}$ and $s$. Based on the analysis of section 3.2, we suggest the following simple local waveform template

$$ \mathcal{F}^{\mathrm{L}}(k_{12}, k_{34}, s) = e^{-\pi\mu} \sin \left[ \mu \left( \arccosh \frac{\gamma k_{12}}{s} - \arccosh \frac{\gamma k_{34}}{s} \right) \right. \\ \left. + \kappa \left( \arccoth \frac{c_s k_{12}}{s} + \arccoth \frac{c_s k_{34}}{s} \right) + \delta_{\mathrm{L}} \right] \quad k_{34}/s \ge c_s^{-1}, \quad (3.75) $$

where $\gamma \equiv c_s/\sqrt{1 - \kappa^2/\mu^2}$ is a dimensionless boost-breaking parameter, and $\delta$ is a constant phase. Similarly, a simple non-local waveform template is given by

<footnote n="10">¹⁰The s-channel dimensionless waveform function only depends on $(k_{12}, k_{34}, s)$ as we consider the simple seed correlator (2.9). The exchange trispectrum corresponding to specific cubic interactions, with a more complicated kinematics dependence, can then be deduced by the application of suitable differential operators, see, e.g., eq. (2.12).</footnote>

<footnote n="11">¹¹Recall that $c_s$ denotes the ratio between the sound speed of inflaton fluctuations and that of the exchanged massive field, and therefore can exceed unity.</footnote>
$$
\mathcal{F}^{\text{NL}}(k_{12}, k_{34}, s) = \begin{cases}
\begin{aligned}
& e^{-\pi(\mu-\kappa)} \cos \left[ \mu \left( \text{arccosh} \frac{\gamma k_{12}}{s} + \text{arccosh} \frac{\gamma k_{34}}{s} \right) \right. \\
& \qquad \left. + \kappa \left( \text{arccoth} \frac{c_s k_{12}}{s} + \text{arccoth} \frac{c_s k_{34}}{s} \right) + \delta_{\text{NL}} \right]
\end{aligned} & k_{34}/s \ge c_s^{-1}, \\
\\
\begin{aligned}
& e^{-\pi(\mu-\kappa)/2} \exp \left[ -\mu \arcsin \frac{\gamma k_{34}}{s} + \kappa \arctan \frac{c_s k_{34} \kappa}{s \mu |\Delta_R|} \right] \\
& \qquad \times \cos \left[ \mu \text{arccosh} \frac{\gamma k_{12}}{s} + \kappa \text{arctanh} \frac{c_s k_{34}}{s} + \delta_{\text{NL}} \right]
\end{aligned} & k_{34}/s \le c_s^{-1}.
\end{cases} \tag{3.76}
$$

with $\Delta_R = \sqrt{c_s^2 k_{34}^2 / s^2 - (1 - \kappa^2 / \mu^2)}$. Let us now comment this template, and show that it encodes all the boost-breaking physics.

* de Sitter limit: in the canonical de Sitter invariant case $\kappa = 0$ and $c_s = 1$, we have $\gamma = 1$ and the edge region ($k_{34}/s \le c_s^{-1}$) disappears due to the triangle inequality. Therefore, the local and non-local signals reduce to

$$
\begin{align}
\mathcal{F}^{\text{L}}(k_{12}, k_{34}, s) &= e^{-\pi\mu} \sin \left[ \mu \left( \operatorname{arccosh} \frac{k_{12}}{s} - \operatorname{arccosh} \frac{k_{34}}{s} \right) + \delta_{\text{L}} \right], \tag{3.77} \\
\mathcal{F}^{\text{NL}}(k_{12}, k_{34}, s) &= e^{-\pi\mu} \cos \left[ \mu \left( \operatorname{arccosh} \frac{k_{12}}{s} + \operatorname{arccosh} \frac{k_{34}}{s} + \delta_{\text{NL}} \right) \right]. \notag
\end{align}
$$

Notice in particular that both the local and the non-local signals have the same Boltzmann suppressed amplitude. In the soft limit $s \to 0$, one recovers the conventional waveforms

$$
\mathcal{F}^{\text{L}} \approx e^{-\pi\mu} \sin[\mu \log(k_{12}/k_{34}) + \delta_{\text{L}}], \quad \mathcal{F}^{\text{NL}} \approx e^{-\pi\mu} \cos[\mu \log(k_{12}k_{34}/s^2) + \delta_{\text{NL}}]. \quad (3.78)
$$

In figure 14, we show both the local and non-local signals in the dimensionless trispectrum waveform function $\mathcal{F}$ in the $s$-channel, varying the mass parameter $\mu$. Remarkably, the refined waveform (3.77) matches well the exact result, even close to regular kinematic configurations $k_{12}/k_{34} \approx 1$ and $k_{12}k_{34}/s^2 \approx 1$, even for $\mu = O(1)$. Meanwhile, the conventional template exhibits an incorrect phase for $\mu \lesssim 3$, and fails to reproduce the amplitude of the local signal as the mass parameter increases. Notice that for the local signal, the conventional template does not accurately reproduce the exact signal in mildly soft kinematic configurations, where the signal-to-noise ratio would be the greater.

* Reduced sound speed limit: the notable feature of a reduced sound speed as the only boost-breaking effect is to modify the amplitude of the non-local signal which reads

$$
\mathcal{F}^{\text{NL}}(k_{12}, k_{34}, s) = \begin{cases} e^{-\pi\mu/2} \cos \left[ \mu \left( \text{arccosh} \frac{c_s k_{12}}{s} + \text{arccosh} \frac{c_s k_{34}}{s} \right) \right] & k_{34}/s \ge c_s^{-1}, \\ e^{-\pi\mu/2 - \mu \arcsin c_s k_{34}/s} \cos \left[ \mu \text{arccosh} \frac{c_s k_{12}}{s} + \delta_{\text{NL}} \right] & k_{34}/s \le c_s^{-1}. \end{cases} \tag{3.79}
$$
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Dimensionless de Sitter invariant trispectrum waveform function $\mathcal{F}(k_{12}, k_{34}, s)$ in the s-channel, as defined in (3.72), as function of the momentum ratio $k_{12}/k_{34}$ for the local signal (left panel), and as function of the momentum ratio $k_{12}k_{34}/s^2$ for the non-local signal (right panel), varying the mass parameter $\mu = 1, 2, 3$. The exact result (black) is expressed in terms of product of hypergeometric functions, see section 2, the refined template (blue) is given by (3.77), and the conventional template (red) is given in (3.78). We have taken $\delta = 0$.</figcaption></figure>

In the ultra-reduced sound speed case $c_s \to 0$ deep in the edge region, this template reduces to the conventional one

$$ \mathcal{F}^{\text{NL}} \approx e^{-\pi\mu/2} \cos[\mu \log(s/c_s k_{12}) + \delta_{\text{NL}}], \qquad (3.80) $$

which describes an amplified cosmological collider signal [11, 26]. Notice that the kinematic-dependent amplitude factor $\exp[-\mu \arcsin(c_s k_{34}/s)]$ interpolates from 1 to $e^{-\pi\mu/2}$ as we scan the edge region $0 < c_s k_{34}/s < 1$. As depicted in figure 15, where we show the dimensionless trispectrum exact, approximate, and conventional waveform, this dependence is crucial to accurately capture the amplitude of the cosmological collider signal. Indeed, we notice that the amplitude of the cosmological collider signal given by the conventional template is inaccurate, especially close to the interior region of the kinematic configurations.

* **Helical template:** further turning on the helical chemical potential yields the full boost-breaking trispectrum templates (3.75) and (3.76). Notably, the non-local signal amplitude is boosted by an enhanced Boltzmann factor $e^{-\pi(\mu-\kappa)/2}$, with again the additional kinematic dependence smoothly interpolating between the conventional cosmological amplitude and that with a reduced sound speed.

**Bispectrum shape template.** The bispectrum can be obtained from the trispectrum by taking the soft limit $k_4 \to 0$. In practice, this amounts to set $k_{34} \to k_3$ and $s \to k_3$ by momentum conservation. The local and non-local signals become indistinguishable, and we
<figure><img src="image_19.png" /><figcaption>Figure 15. Dimensionless trispectrum waveform function $\mathcal{F}_s(k_{12}, k_{34}, s)$ in the s-channel, as defined in (3.72), as function of the momentum ratio $k_{12}/k_{34}$ in the case of a reduced sound speed $c_s = 0.1$ in the edge region (left panel), and further turning on the chemical potential $\kappa \neq 0$ (right panel), varying the mass parameter $\mu = 1, 2, 3$. The solid lines correspond to $k_{34}/s = 1$ and the dashed line corresponds to $k_{34}/s = 5$, closer to the interior region. We refer to figure 14 for details about the exact, approximate, and conventional templates. We have taken $\delta = 3\pi/4$. The shaded grey region corresponds to the diamond region in red shown in figure 10.</figcaption></figure>

refer to the resulting signal collectively as the cosmological collider signal. Following standard conventions, we define the dimensionless bispectrum shape function $S$ such that

$$ \langle \zeta_{k_1} \zeta_{k_2} \zeta_{k_3} \rangle = (2\pi)^7 \delta^{(3)}(\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3) \frac{\mathcal{P}_{\zeta}^2}{(k_1 k_2 k_3)^2} S(k_1, k_2, k_3). \quad (3.81) $$

Based on the previous dimensionless trispectrum, we propose the following boost-breaking bispectrum shape template

$$ S^{\text{collider}}(k_1, k_2, k_3) = \mathcal{B}(\{k_i\}) P_S(\hat{\mathbf{k}}_1 \cdot \hat{\mathbf{k}}_3) e^{-\pi\mu/2 - \mu \arcsin c_s} \\ \times \cos \left[ \mu \operatorname{arccosh} \left( \frac{c_s k_{12}}{k_3} \right) + \delta \right] + \text{perms.}, \quad (3.82) $$

where $\delta$ is $\mu$-dependent constant phase, and $\mathcal{B}(\{k_i\})$ is a form factor that depends on the considered cubic interaction, and in general has a polynomial dependence on the mass parameter $\mu$. As an example, the form factor from the $\dot{\varphi}^2\sigma - \dot{\varphi}\sigma$ scalar exchange reads $\mathcal{B}(\{k_i\}) = \mu^{3/2}k_1k_2/(k_{12})^2[k_3/(k_{12})]^{1/2}$. Other precise momentum dependences can be found in [6]. Note that in the ultra-soft limit $k_3 \to 0$, one recovers the standard functional dependence $S^{\text{collider}} \sim \cos[\mu \log(k_{12}/k_3) + \delta]$. The template (3.82) is more accurate in the intermediate parameter range $\mu = \mathcal{O}(1)$ and $c_s \lesssim 1$, and also matches well the cosmological collider signal in mildly soft kinematic configurations. Notice that the bispectrum requires a linear mixing between the external and the exchanged massive fields, such as $\partial_{i_1...i_S} \varphi \sigma_{i_1...i_S}$. Therefore, only the longitudinal mode of the spinning particle can contribute to the bispectrum, while other
components with $\lambda \neq 0$ are projected out by the linear mixing. This leads to the universal angular dependence $P_S(\hat{k}_1 \cdot \hat{k}_3)$ written in (3.82) and no chemical potential. This refined cosmological collider template is both straightforward to implement and practically valuable. The fact that it consists solely of elementary functions significantly reduces computational costs. At the same time, it serves as an excellent approximation to the exact result in mildly soft kinematic regimes, where the signal-to-noise ratio is highest. Finally, it provides a smooth interpolation between the de Sitter invariant and strongly boost-breaking limits.

# 4 Conclusions

In this paper, we have computed the most general tree-level boost-breaking seed correlator and have classified signals of cosmological collider physics based on our newly developed saddle-point method to evaluate bulk time integrals. The resulting analytical shapes offer valuable templates for cosmological collider studies. These template waveforms are particularly relevant for phenomenologically motivated scenarios with potential observational signatures in current and upcoming cosmological surveys.

More specifically, we have focused on cosmological correlators sourced by the exchange of (spinning) particles and external inflatons with reduced speed of sound. This technically extends the kinematically allowed configurations where standard series solutions fail to converge. To compute the correlators, we employ both the well-established bootstrap approach and a newly developed spectral method, each offering complementary advantages in terms of convergence and computational efficiency. We then provide a detailed and handy recipe to extract and classify various cosmological collider signals based on expanding the internal mode using a WKB approximation and evaluating the time integrals with the saddle-point method. Notably, the saddle points in the complex time plane map directly onto the dynamics of massive fields, providing clear physical understanding of various bulk processes. We then showcase this method by studying in detail de Sitter-invariant and boost-breaking scenarios. The study on the motion of the saddle point in the complex plane driven by external kinematics reveals a deep connection between the amplitude of boundary kinematical discontinuities and the residue of bulk time singularities, giving a natural explanation to the fractional Boltzmann suppression of cosmological collider signals. By virtue of the saddle-point analysis, we also find novel new signals — e.g. oscillating linearly in the momentum ratio with a frequency set by the product of the speed of sound and the chemical potential appearing in the edge regions—, as well as a refined waveform template for the standard cosmological collider signals. Importantly, after all boost-breaking effects are captured at tree-level, the resulting signals can be sizeable and potentially lie within reach of current or near-future cosmological observations.

Eventually, our work opens several interesting avenues for future investigation:

* **Linking kinematic and time singularities.** In section 3, we noticed an intriguing connection between the amplitudes of the boundary kinematic singularities and poles in the bulk time. In all cases encountered so far, the discontinuity along the branch cut of kinematic variables are dictated by the exponential of the residue of the dispersion relation at time-zero or infinity. We believe such a connection is no mere coincidence,
but reflects a deep correspondence between the boundary and the bulk, which certainly
deserves further exploration in the future.

* **Towards saddle points in loops.** The saddle-point analysis for cosmological collider signals we adopted in this paper is exclusively focused on tree diagrams. At loop level, one needs to incorporate momentum integrals in addition to time integrals. Given that the cutting rule still applies in loop diagrams [51, 81], it is conceivable that a similar saddle-point analysis can be performed there, yielding simplified shape templates for cosmological collider signals at loop level. We leave this interesting possibility for future work.

* **Breaking scale invariance.** In this work, we computed the most general boost-breaking correlator using different of analytical and approximation methods. An important and complementary direction for future research involves breaking of scale invariance (also known as the dS dilation), which is commonly violated in many inflationary models with features. Such breaking of scale invariance often leads to enhanced cosmological collider signals [23, 25, 28, 39, 40]. It is therefore of significant interest to develop a general analytical framework for studying these symmetry-breaking scenarios, as well as to extend our approximation methods to gain deeper physical insights into such cases.

* **Searching with refined templates.** Our approximation scheme produces many refined templates for cosmological collider signals that are ready to be searched for in present and upcoming data. In particular, theories with broken de Sitter boosts naturally allow for a large amplitude of the trispectrum, over which we now have precise and efficient analytical control. It is thus interesting to put our templates to the test.

## Acknowledgments

We thank Shuntaro Aoki, Xingang Chen, Thomas Colas, Sadra Jazayeri, Fumiya Sano, David Stefanyszyn, Yi Wang, Dong-Gang Wang, Zhong-Zhi Xianyu, Masahide Yamaguchi, and Bowei Zhang for helpful discussions. Z.Q. is supported by NSFC under Grant No. 12275146, the National Key R&D Program of China (2021YFC2203100), and the Dushi Program of Tsinghua University. During the course of this work, S.R-P. and D.W. were supported by the European Research Council under the European Union's Horizon 2020 research and innovation programme (grant agreement No 758792, Starting Grant project GEODESI). X.T. is supported by STFC consolidated grants ST/T000694/1 and ST/X000664/1. Y.Z. is supported by the IBS under the project code, IBS-R018-D3. This article is distributed under the Creative Commons Attribution International Licence (CC-BY 4.0).

## A Mathematical interlude

Classifying non-analyticities of cosmological correlators employs tools from complex analysis. To ensure this paper remains self-contained and accessible, this appendix offers a pedagogical and elementary introduction to the theory of Riemann surfaces and the Saddle point approximation. Emphasis is placed on clarity and intuition rather than formal mathematical rigour.
## A.1 Riemann surfaces from multi-valued functions

The study of differentiable functions of a single complex variable naturally unveils the concept of multi-valued functions. This, in turn, introduces an elegant class of mathematical objects known as Riemann surfaces, which beautifully connects analysis to geometry.<fnref n="12" /> The intuitive idea behind Riemann surfaces is to take multiple branches of a complex multi-valued function and glue these together to obtain a continuous single-valued function. The main question we will try answer is: “What do we really mean by multi-valued functions?”.

**Complex logarithm.** Let us take the canonical example of the complex logarithm function to introduce the concept of multi-valued functions and the corresponding Riemann surface. We try to extend the real logarithm $x \in \mathbb{R}_+^* \mapsto \log(x)$ to the largest possible open, connected subset of the complex plane $\mathbb{C}$. The most natural way to proceed is to come back to the defining property of the real logarithmic function: it is the primitive of $x \in \mathbb{R}^* \mapsto 1/x$ that vanishes at $x = 1$. We therefore define

$$ \log(z) = \int_{C(z)} \frac{dz'}{z'}, \quad \text{for } z \in \mathbb{C}^*, \qquad (A.1) $$

where $C(z)$ is a contour linking 1 to $z$ in the complex plane. Eq. (A.1) defines a function if every choice of the contour $C(z)$ leads to the same value for $\log(z)$. Let us consider two contours $C_1(z)$ and $C_2(z)$ linking 1 to $z$, but passing on either side of the origin, as depicted on the right.

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

If we denote $C$ the closed contour composed of $C_1$ and $-C_2$, then using the residue theorem, we obtain $\int_{C(z)} \frac{dz'}{z'} = 2i\pi$. This yields

$$ \int_{C_1(z)} \frac{dz'}{z'} = \int_{C_2(z)} \frac{dz'}{z'} + 2i\pi. \qquad (A.2) $$

It is therefore not possible to define the complex logarithm on the entire punctured complex plane $\mathbb{C}\setminus\{0\}$. To avoid such inconvenience, we restrict the argument of $z$, that we denote $\arg(z)$, to $\arg(z) \in (\alpha, \alpha + 2\pi)$ with real $\alpha$. The line $\{z | \arg(z) = \alpha\}$ is a branch cut, with the origin $z = 0$ being the branch point, i.e. the common point to all branch cuts. Choosing $\alpha = -\pi$ so that the branch cut is the negative real axis $\mathbb{R}_-$ defines the principal branch of the complex logarithm. However, note that any other value for $\alpha$, defining a line $C$ linking the origin to infinity, would give a valid branch cut and a new complex logarithm can be defined on $C\setminus C$. A natural consequence is if we consider two complex logarithms $\ell(z)$ and $\ell'(z)$ defined on $C\setminus C$ and $C\setminus C'$, respectively, the difference $\ell(z) - \ell'(z)$ is a multiple of $2i\pi$ for all $z \in C\setminus\{C \cup C'\}$.

<footnote n="12">¹²For a general introduction to Riemann surfaces, see e.g. [90].</footnote>
Importantly, a different point of view removes the need to introduce a branch cut to define the complex logarithm. In fact, when physicists' point of view is that of a multi-valued function, the mathematicians' one is that of a function defined on a Riemann surface.

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

Introducing a Riemann surface corresponds to replacing the punctured complex plane $C\setminus\{0\}$, on which we tried without succeeding to define a complex holomorphic logarithm, by a larger connected domain. Intuitively, starting from the principal branch, rotating counter-clockwise once around the origin would pick a factor $2i\pi$ to the logarithm to end up on the next branch, and so on. We therefore can construct the following surface: take an infinite number of complex planes that we stack on top of each other. Each of these planes are cut on a branch cut, say the negative real axis $\mathbb{R}_-$, and the upper branch cut is glued to the lower one. The final result is a connected surface — a helicoid — resembling a spiralling parking garage with infinitely many levels, as shown on the right: it is the *Riemann surface* of the complex logarithm where the z-axis is the imaginary part of $\log(z)$. A point on this Riemann surface $\mathfrak{R}$ can be thought of as a pair $(z, n)$, with $n \in \mathbb{Z}$ corresponding to a Riemann sheet for which the complex logarithm would give $\log(r) + i(\theta + 2\pi n)$, where $z = re^{i\theta}$ with $r \in \mathbb{R}_+^*$ and $\theta = \arg(z)$. As such, the Riemann surface can be embedded in $C \times \mathbb{Z}$. On $\mathfrak{R}$, the complex logarithm can be defined as follows: (i) we choose a contour $C$ linking $(1, 0) \in \mathfrak{R}$ to the point $(z, n) \in \mathfrak{R}$, (ii) this contour can be projected vertically on the principal branch $n = 0$, resulting in a contour joining 1 to $z$, and (iii) we define $\log(z) = \int_C dz'/z'$. This way, the function does not depend on the contour any more as whenever we circle around the origin, we move to another branch and the logarithm picks up a factor $\pm 2i\pi$.

**Complex square root.** Another classic example of a function that presents challenges when extending its domain from the positive real axis to the complex plane is the square root. Indeed, defining $z = re^{i\theta}$ with $r \in \mathbb{R}_+^*$ and $\theta = \arg(z)$, a naive definition of the square root would simply be $\sqrt{z} = \sqrt{r}e^{i\theta/2}$. However, the argument $\theta/2$ is only defined up to a factor $\pi$, resulting in a sign ambiguity on the sign of $\sqrt{z}$. This function is multi-valued: for a single value of $z$, there are two different choices for $\sqrt{z}$. We define these two choices by

$$f_>(z) = \sqrt{r}e^{i\theta/2}, \quad f_<(z) = \sqrt{r}e^{i(\theta/2+\pi)}, \qquad (A.3)$$

where $f_>$ and $f_<$ are the positive and negative branches of the square root, respectively. As we observed earlier with the complex logarithm, it is necessary to define a branch cut. From now on, we will adopt the principal branch cut, positioned along the negative real axis $\mathbb{R}_-$, i.e. $\theta = \arg(z) \in (-\pi, \pi)$. Let us now look at the behaviour of each of the two branches
$f_>(z)$ and $f_<(z)$ as $z$ approaches the branch cut from above and from below on a circle centred at the origin with a radius $r > 0$. Taking the limit $\theta \to +\pi$, we obtain $f_>(z) = +i\sqrt{r}$ and $f_<(z) = -i\sqrt{r}$. Similarly, as $\theta \to -\pi$, we have $f_>(z) = -i\sqrt{r}$ and $f_<(z) = +i\sqrt{r}$. We represent these limits below for the positive and negative branches.

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

To define a function on the punctured complex plane $C\setminus\{0\}$ that remains continuous along the branch cut, we observe that the sign is reversed after a single rotation around the origin but returns to its original value after two full rotations. We can e.g. approach the branch cut from above using the positive branch $f_>(z)$ and exit from below using the negative branch $f_<(z)$. The corresponding Riemann surface for the complex square root is depicted on the right: when the square root takes a specific value on one Riemann sheet, it assumes the opposite value on the other sheet. We end up with a self-intersecting surface.

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

**Constructing Riemann surfaces.** Whenever a function is multi-valued, it can be promoted to a single-valued function albeit extending the domain of validity to the corresponding Riemann surface. For pedagogical purposes, let us now explicitly construct the Riemann surface corresponding to the function describing the time-dependent frequency of a massive field scalar field in de Sitter: (dimensionless) frequencies $\omega$ such as

$$ \omega^2(z) = 1 + \frac{1}{z^2}. \qquad (A.4) $$

The function $\omega(z)$ has a pole at the origin $z = 0$ and two branch points at $z = \pm i$. Taking the square root principal branch, we place the branch cut on the imaginary axis between $+i$ and $-i$. When trying to extend the domain of validity of the frequency $\omega(z)$ to the largest possible domain of the complex plane of $z$, we note that this function has two branches corresponding to the positive and negative branches of the square root, that we denote

$$ \omega_>(z) = +\sqrt{1 + \frac{1}{z^2}}, \quad \omega_<(z) = -\sqrt{1 + \frac{1}{z^2}}, \qquad (A.5) $$

that are defined on two distinct Riemann sheets. We now probe the behaviour of both branches on each sides of both branch cuts. On a small circle with infinitesimal radius $\epsilon > 0$ centred at
the origin, we have $\omega_>(\epsilon e^{i\theta}) \approx e^{-i\theta}/\epsilon \to \mp i/\epsilon$ as $\theta \to \pm\pi/2$. Similarly for the negative branch, we obtain $\omega_<(\epsilon e^{i\theta}) \approx e^{-i(\theta-\pi)}/\epsilon \to \pm i/\epsilon$ as $\theta \to \pm\pi/2$, which can be visualised as follows:

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

The Riemann surface can be constructed by noticing that we can continuously enter the upper portion of the branch cut from the left using the positive branch to exit the lower portion of the branch cut using the negative branch cut. However, when crossing the imaginary axis away from the branch cut, we remain on the same Riemann sheet. We obtain the Riemann surface used in the main text.

## A.2 Saddle point approximation of complex integrals

Local behaviours of integrals in certain limits can completely determine their asymptotic forms without introducing free parameters. This is in sharp contrast with a local study of differential equations around singular points where solutions are always found up to arbitrary integration constants. Among these local analysis, the saddle point approximation — also called the method of steepest descent — is a method to find approximate solutions to complex integrals where one deforms the integration contour to pass a saddle (stationary) point in the direction of steepest descent or stationary phase. We invite the reader to consult e.g. [91, 92] for more details.

**Preliminary example.** The main idea of the saddle point approximation is to deform the integration contour in the complex plane so that the integrand has a *constant phase* on the new contour. Before formally exposing the method, we consider a preliminary example to motivate why it is helpful to shift complex contours to perform an asymptotic analysis. Let us consider the following integral:

$$I(\mu) = \int_{0}^{1} dt \log(t) e^{i\mu t}. \qquad (A.6)$$

As $\mu \gg 1$, the asymptotic expansion of $I(\mu)$ is non-analytic, which complicates finding an approximate solution. We therefore deform the integration interval in blue to the complex
contour $C = C_1 + C_2 + C_3$ as depicted in red below:

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

By Cauchy's theorem, we have $I(\mu) = \int_{C_1+C_2+C_3} dz \log(z) e^{i\mu z}$. As $t_* \to \infty$, the integral over $C_2$ vanishes exponentially. The integral (A.6) becomes

$$I(\mu) = i \int_{0}^{\infty} dz \log(iz) e^{-\mu z} - i \int_{0}^{\infty} dz \log(1 + iz) e^{i\mu(1+iz)}. \quad (A.7)$$

The first integral can be performed exactly, yielding

$$i \int_{0}^{\infty} dz \log(iz) e^{-\mu z} = -i \frac{\log(\mu)}{\mu} - \frac{i\gamma + \pi/2}{\mu}, \quad (A.8)$$

where we have used the principal branch of the complex logarithm $\log(iz) = \log(z) + i\pi/2$ and $\int_0^\infty \log(t)e^{-t}dt = -\gamma$ with $\gamma = 0.5772...$ the Euler's constant. The second integral can also be performed by Taylor expanding the logarithm and integrating each term in the series

$$i \int_{0}^{\infty} dz \log(1 + iz) e^{i\mu(1+iz)} = i e^{i\mu} \sum_{n=1}^{\infty} \frac{(-i)^n (n-1)!}{\mu^{n+1}}. \quad (A.9)$$

Importantly, notice that along the contours $C_1$ and $C_3$, the phase of $e^{i\mu t}$ is constant, with frequency 0 and 1, respectively. These contours are also the steepest-descent paths as $|e^{i\mu t}|$ decreases most rapidly (in fact exponentially) as $t$ ranges from 0 and 1 to $\infty$. This is no coincidence as we now review.

**Steepest-descent complex contours passing through saddle points.** The method of steepest descent enables one to approximate complex integrals of the form

$$I(\mu) = \int_C dz g(z) e^{\mu f(z)}, \quad (A.10)$$

in the limit $\mu \gg 1$. The contour $C$ can be deformed into a new contour $C'$ without changing the integral provided no poles nor branch cuts are crossed. If one chooses a path that satisfies the following conditions: (i) the phase $\operatorname{Im}[f(z)] \equiv \theta$ is constant along $C'$ and (ii) the contour $C'$ passes through one or more saddle points $z_\bullet$ such that $f'(z_\bullet) = 0$, then the integral reads

$$I(\mu) = e^{i\mu \operatorname{Im}[f(z)]} \int_{C'} dz g(z) e^{\mu \operatorname{Re}[f(z)]} \approx e^{i\mu\theta} \sqrt{\frac{2\pi}{\mu \frac{d^2}{dz^2} \operatorname{Re}[f]|_{z_\bullet}}} g(z_\bullet) e^{\mu \operatorname{Re}[f(z_\bullet)]}, \quad (A.11)$$

where the remaining integral is approximated using Laplace's method, which follows from Taylor expanding the function $\operatorname{Re}[f]$ around the stationary point $z = z_\bullet$ and finishing the
Gaussian integral by extending the limits of integration by $\pm\infty$.<fnref n="13" /> Extensions to arbitrary precision can be systematically found. In the presence of multiple saddles, one needs to sum their contributions. However, it is in general expected that the saddle points that are located the closest to the original integration contour contribute the most. Note that if the function $f$ is analytic, constant phase contours are equivalent to steepest-descent contours. Indeed, setting $z = x + iy$, using the Cauchy-Riemann equations, we obtain

$$ \frac{\partial \operatorname{Re}[f]}{\partial x} \frac{\partial \operatorname{Im}[f]}{\partial x} + \frac{\partial \operatorname{Re}[f]}{\partial y} \frac{\partial \operatorname{Im}[f]}{\partial y} = 0 \quad \Leftrightarrow \quad \nabla \operatorname{Re}[f] \cdot \nabla \operatorname{Im}[f] = 0, \qquad (A.12) $$

which shows that $\operatorname{Im}[f]$ is constant on contours whose tangents are parallel to $\nabla \operatorname{Re}[f]$. Two distinct steepest curves can only intersect at saddle points.

**Illustrative example.** We now apply this method to a concrete example. Let us derive the asymptotic expansion of the integral

$$ I(\mu) = \int_C \frac{ds}{2i\pi} e^s s^{-\mu}, \qquad (A.13) $$

in the limit $\mu \gg 1$, where $C$ is the Hankel contour that begins at $-\infty - i\epsilon$, encircles the branch cut lying on the real negative axis and ends up at $-\infty + i\epsilon$ with $\epsilon > 0$. We first make the change of variables $s = \mu z$ to convert the movable saddle point that depends on $\mu$ to a fixed one, yielding

$$ I(\mu) = \frac{1}{\mu^{\mu-1}} \int_C \frac{dz}{2i\pi} e^{\mu(z-\log z)}. \qquad (A.14) $$

Defining $f(z) = z - \log z$, solving $f'(z) = 0$ gives a single saddle point at $z_\bullet = 1$. To study the structure of this saddle, we set $z = x + iy$ so that $f(z) = x - \log \sqrt{x^2 + y^2} + i(y - \arctan y/x)$. Therefore, paths of constant phase emerging from $z_\bullet = 1$ (as $f(z_\bullet) = 1$) satisfy

$$ y = \arctan y/x. \qquad (A.15) $$

There are two solutions to this equation: $y = 0$ and $x = y/\tan(y)$. Now we can deform the original contour $C$ in blue to $C'$ represented in red that passes by the saddle.<fnref n="14" />

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

<footnote n="13">¹³This creates only a small error as the exponential decays very fast away from the saddle point.</footnote>

<footnote n="14">¹⁴It is not possible to continuously deform the original contour to the steepest-descent path spanning the positive real axis.</footnote>
The constant-phase (or steepest-descent) curve $C' = \{(t/\tan(t), t)|t \in \mathbb{R}\}$ is parallel to $\operatorname{Im}[z] = y$ around $z_\bullet = 1$ and crosses the imaginary axis at $\pm i\pi/2$. Along this path, $|e^{\mu(z-\log z)}|$ reaches its maximum at the saddle and decays exponentially as $z \to -\infty \pm i\pi$. In the case of no saddle point, $|e^{\mu(z-\log z)}|$ would decrease in one direction along $C'$ and increase in the other direction, and the integral would not be convergent. Note that if instead the steepest-descent path is finite and does not pass through a saddle point, then the maximum of the integrand must occur at an endpoint of the contour. To finish the integral, we examine the neighbourhood of $z_\bullet$ along $C'$ by setting $z = 1 + it$ with $t \ll 1$, obtaining

$$I(\mu) \approx \frac{e^{\mu}}{\mu^{\mu-1}} \int_{-\epsilon}^{+\epsilon} \frac{dt}{2\pi} e^{-t^2/2} \approx \frac{e^{\mu}}{\mu^{\mu-1}} \frac{1}{\sqrt{2\pi\mu}}, \qquad (A.16)$$

where we have let $\epsilon \to \infty$. Note that it is important to stay on the constant-phase contour as naively applying (A.11) would have given an overall wrong phase $e^{-i\pi/2}$.

**Application to massive field in de Sitter.** We now treat in detail the simplest case of a massive scalar field in de Sitter with zero chemical potential $\kappa = 0$ and unity sound speed $c_s = 1$, studied in section 3.2.1. After applying the leading-order WKB approximation to the massive mode function and discarding the negative-frequency integral as it is Boltzmann suppressed, we are left with the following integral

$$F_L^{(\alpha)} = \int_{-\infty+}^{0} \frac{d\tau}{\sqrt{2\omega_s(\tau)}} e^{i(k_{12}\tau - \int_{\tau_i}^{\tau} \omega_s(\tau')d\tau')}, \qquad (A.17)$$

where $\tau_i$ is an arbitrary time in the infinite past, and the effective frequency $\omega_s(\tau)$ is defined by

$$\omega_s^2(\tau) = s^2 + \frac{\mu^2}{\tau^2}. \qquad (A.18)$$

The saddle points are found to $\tau_\bullet = \pm\mu/\sqrt{k_{12}^2 - s^2}$ which lie on the real axis in the complex $\tau$ plane as $k_{12} > s$ by the triangle inequality. The closest saddle point to the original integration contour is the one located on the negative real axis, so we expect its contribution to dominate the integral. Let us now evaluate the Dingle's singulant function

$$F_s(\tau) = -2i \int_{\tau_c}^{\tau} \omega_s(\tau')d\tau', \qquad (A.19)$$

with $\tau_c = +i\mu/s$ the critical point (with negative imaginary part) that makes the effective frequency vanish. Selecting the positive branch of the square root, we obtain

$$F_s(\tau) = -2i \left[ \tau \sqrt{s^2 + \mu^2/\tau^2} + \mu \operatorname{arctanh} \left( \frac{-\mu}{\tau \sqrt{s^2 + \mu^2/\tau^2}} \right) \right]_{-i\mu/s}^{\tau} \qquad (A.20) \\ = -2i \left( \tau \sqrt{s^2 + \mu^2/\tau^2} + \mu \operatorname{arctanh} \left( \frac{-\mu}{\tau \sqrt{s^2 + \mu^2/\tau^2}} \right) \right) + \pi\mu.$$

After splitting the integral as $\int_{\tau_c}^{\tau} = \int_{\tau_i}^{\tau} - \int_{\tau_i}^{\tau_c}$, the exponent in the integral (A.17) therefore can be written in terms of $F(\tau)$ as

$$f(\tau) = ik_{12}\tau + \frac{1}{2}F_s(\tau) - \pi\mu/2, \qquad (A.21)$$
where we have discarded the contribution coming from $\tau_i$ as it is arbitrary. Let us now study the steepest-descent curve passing through the saddle $\tau_\bullet$. To do that, we Taylor expand the exponent in the neighbourhood of the saddle $\tau_\bullet$

$$f(\tau_\bullet + \epsilon) \approx f(\tau_\bullet) - i\alpha\epsilon^2, \qquad (A.22)$$

with $|\epsilon| \ll 1$, where $f(\tau_\bullet) = -i\mu \operatorname{arctanh} \sqrt{k_{12}^2 - s^2}/k_{12} \equiv i\beta$ and $\alpha \equiv (k_{12}^2 - s^2)^{3/2}/2k_{12}\mu$. We set $\epsilon = x + iy$ and obtain

$$f(\tau_\bullet + \epsilon) \approx 2\alpha xy + i[\beta - \alpha(x^2 - y^2)]. \qquad (A.23)$$

Since $f(\tau_\bullet)$ is purely imaginary with a negative imaginary component, the constant-phase curve crossing the saddle must satisfy the following equation

$$\beta - \alpha(x^2 - y^2) = -\lambda, \qquad (A.24)$$

where $\lambda > 0$ is a free parameter that labels the different constant-phase curves. For $\tau_\bullet < 0$, the solution to this equation is $x = -\sqrt{y^2 + (\beta + \lambda)/\alpha}$, and the constant-phase curves in the neighbourhood of the saddle are

$$\gamma_{\lambda}(t) = \{(-\sqrt{t^2 + (\beta + \lambda)/\alpha}, t) | t \in \mathbb{R}\}. \qquad (A.25)$$

The affine parameter where the curve crosses the saddle is $t_\bullet = 0$ and the label selecting this particular curve is found to be $\lambda = \mu\sqrt{k_{12}^2 - s^2}/2k_{12} + \mu \operatorname{arctanh} \sqrt{k_{12}^2 - s^2}/k_{12}$. The tangent vector at the saddle point is found to be $T = (0, 1)$. Therefore, the curve is parallel to the imaginary axis. We illustrate the complex $\tau$ plane below, where the two saddle points are depicted in black, the critical points in red along with the branch cut joining them, and the constant-phase curve passing by the saddle in blue (whose expression is only found in the neighbourhood of $\tau_\bullet$). Additional constant-phase curves are shown in gray.

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

To expand the exponent (A.21) around the saddle point along the steepest-descent path, we set $\tau = \tau_\bullet - it$ with $t \ll 1$ as the original contour starts at $-\infty^+$. We obtain

$$F_L^{(\alpha)} = \frac{-ie^{-i\mu} \operatorname{arctanh} \frac{\sqrt{k_{12}^2 - s^2}}{k_{12}}}{\sqrt{2\omega_s(\tau_\bullet)}} \int_{-\epsilon}^{+\epsilon} dt e^{i \frac{(k_{12}^2 - s^2)^{3/2}}{2k_{12}\mu} t^2} \qquad (A.26) \\ \approx \frac{\sqrt{\pi\mu}}{(k_{12}^2 - s^2)^{3/4}} e^{-i\pi/4} \exp \left[ -i\mu \operatorname{arccosh} \frac{k_{12}}{s} \right],$$
where we have finished the Gaussian integral after setting $\epsilon \to \infty$. Note that integrating along the constant-phase curve is essential to obtain the correct phase.

**Data Availability Statement.** This article has associated data in a data repository. Available at: https://github.com/deniswerth/Boost-Breaking-Cosmological-Correlators

## References

[1] P.D. Meerburg et al., *Primordial Non-Gaussianity*, Bull. Am. Astron. Soc. **51** (2019) 107 [arXiv:1903.04409] [INSPIRE].

[2] A. Achúcarro et al., *Inflation: Theory and Observations*, arXiv:2203.08128 [INSPIRE].

[3] X. Chen and Y. Wang, *Large non-Gaussianities with Intermediate Shapes from Quasi-Single Field Inflation*, Phys. Rev. D **81** (2010) 063511 [arXiv:0909.0496] [INSPIRE].

[4] X. Chen and Y. Wang, *Quasi-Single Field Inflation and Non-Gaussianities*, JCAP **04** (2010) 027 [arXiv:0911.3380] [INSPIRE].

[5] D. Baumann and D. Green, *Signatures of Supersymmetry from the Early Universe*, Phys. Rev. D **85** (2012) 103520 [arXiv:1109.0292] [INSPIRE].

[6] T. Noumi, M. Yamaguchi and D. Yokoyama, *Effective field theory approach to quasi-single field inflation and effects of heavy fields*, JHEP **06** (2013) 051 [arXiv:1211.1624] [INSPIRE].

[7] N. Arkani-Hamed and J. Maldacena, *Cosmological Collider Physics*, arXiv:1503.08043 [INSPIRE].

[8] J.-O. Gong, S. Pi and M. Sasaki, *Equilateral non-Gaussianity from heavy fields*, JCAP **11** (2013) 043 [arXiv:1306.3691] [INSPIRE].

[9] X. Chen, M.H. Namjoo and Y. Wang, *Quantum Primordial Standard Clocks*, JCAP **02** (2016) 013 [arXiv:1509.03930] [INSPIRE].

[10] X. Chen, Y. Wang and Z.-Z. Xianyu, *Standard Model Background of the Cosmological Collider*, Phys. Rev. Lett. **118** (2017) 261302 [arXiv:1610.06597] [INSPIRE].

[11] H. Lee, D. Baumann and G.L. Pimentel, *Non-Gaussianity as a Particle Detector*, JHEP **12** (2016) 040 [arXiv:1607.03735] [INSPIRE].

[12] H. An, M. McAneny, A.K. Ridgway and M.B. Wise, *Quasi Single Field Inflation in the non-perturbative regime*, JHEP **06** (2018) 105 [arXiv:1706.09971] [INSPIRE].

[13] S. Lu, Y. Wang and Z.-Z. Xianyu, *A Cosmological Higgs Collider*, JHEP **02** (2020) 011 [arXiv:1907.07390] [INSPIRE].

[14] A. Hook, J. Huang and D. Racco, *Searches for other vacua. Part II. A new Higgstory at the cosmological collider*, JHEP **01** (2020) 105 [arXiv:1907.10624] [INSPIRE].

[15] T. Liu, X. Tong, Y. Wang and Z.-Z. Xianyu, *Probing P and CP Violations on the Cosmological Collider*, JHEP **04** (2020) 189 [arXiv:1909.01819] [INSPIRE].

[16] S. Kumar and R. Sundrum, *Cosmological Collider Physics and the Curvaton*, JHEP **04** (2020) 077 [arXiv:1908.11378] [INSPIRE].

[17] L.-T. Wang and Z.-Z. Xianyu, *In Search of Large Signals at the Cosmological Collider*, JHEP **02** (2020) 044 [arXiv:1910.12876] [INSPIRE].

[18] L.-T. Wang and Z.-Z. Xianyu, *Gauge Boson Signals at the Cosmological Collider*, JHEP **11** (2020) 082 [arXiv:2004.02887] [INSPIRE].
[19] Y. Wang and Y. Zhu, *Cosmological Collider Signatures of Massive Vectors from Non-Gaussian Gravitational Waves*, JCAP 04 (2020) 049 [arXiv:2001.03879] [INSPIRE].

[20] L. Pinol, S. Aoki, S. Renaux-Petel and M. Yamaguchi, *Inflationary flavor oscillations and the cosmic spectroscopy*, Phys. Rev. D 107 (2023) L021301 [arXiv:2112.05710] [INSPIRE].

[21] Y. Cui and Z.-Z. Xianyu, *Probing Leptogenesis with the Cosmological Collider*, Phys. Rev. Lett. 129 (2022) 111301 [arXiv:2112.10793] [INSPIRE].

[22] M. Reece, L.-T. Wang and Z.-Z. Xianyu, *Large-field inflation and the cosmological collider*, Phys. Rev. D 107 (2023) L101304 [arXiv:2204.11869] [INSPIRE].

[23] X. Chen, R. Ebadi and S. Kumar, *Classical cosmological collider physics and primordial features*, JCAP 08 (2022) 083 [arXiv:2205.01107] [INSPIRE].

[24] Z. Qin and Z.-Z. Xianyu, *Phase information in cosmological collider signals*, JHEP 10 (2022) 192 [arXiv:2205.01692] [INSPIRE].

[25] D. Werth, L. Pinol and S. Renaux-Petel, *Cosmological Flow of Primordial Correlators*, Phys. Rev. Lett. 133 (2024) 141002 [arXiv:2302.00655] [INSPIRE].

[26] S. Jazayeri, S. Renaux-Petel and D. Werth, *Shapes of the cosmological low-speed collider*, JCAP 12 (2023) 035 [arXiv:2307.01751] [INSPIRE].

[27] S. Jazayeri et al., *Parity violation from emergent nonlocality during inflation*, Phys. Rev. D 108 (2023) 123523 [arXiv:2308.11315] [INSPIRE].

[28] L. Pinol, S. Renaux-Petel and D. Werth, *The cosmological flow: a systematic approach to primordial correlators*, JCAP 02 (2025) 019 [arXiv:2312.06559] [INSPIRE].

[29] Y. Ema and S. Verner, *Cosmological collider signatures of Higgs-$R^2$ inflation*, JCAP 04 (2024) 039 [arXiv:2309.10841] [INSPIRE].

[30] X. Tong, Y. Wang, C. Zhang and Y. Zhu, *BCS in the sky: signatures of inflationary fermion condensation*, JCAP 04 (2024) 022 [arXiv:2304.09428] [INSPIRE].

[31] P. Chakraborty and J. Stout, *Compact scalars at the cosmological collider*, JHEP 03 (2024) 149 [arXiv:2311.09219] [INSPIRE].

[32] D. Werth, L. Pinol and S. Renaux-Petel, *CosmoFlow: Python Package for Cosmological Correlators*, Class. Quant. Grav. 41 (2024) 175015 [arXiv:2402.03693] [INSPIRE].

[33] S. Aoki, A. Ghoshal and A. Strumia, *Cosmological collider non-Gaussianity from multiple scalars and $R^2$ gravity*, JHEP 11 (2024) 009 [arXiv:2408.07069] [INSPIRE].

[34] Y.-P. Wu, *The cosmological collider in $R^2$ inflation*, JCAP 07 (2024) 010 [arXiv:2404.05031] [INSPIRE].

[35] N. Craig, S. Kumar and A. McCune, *An effective cosmological collider*, JHEP 07 (2024) 108 [arXiv:2401.10976] [INSPIRE].

[36] C. McCulloch, E. Pajer and X. Tong, *A cosmological tachyon collider: enhancing the long-short scale coupling*, JHEP 05 (2024) 262 [arXiv:2401.11009] [INSPIRE].

[37] G. Cabass et al., *BOSS constraints on massive particles during inflation: The cosmological collider in action*, Phys. Rev. D 111 (2025) 063510 [arXiv:2404.01894] [INSPIRE].

[38] W. Sohn, D.-G. Wang, J.R. Fergusson and E.P.S. Shellard, *Searching for cosmological collider in the Planck CMB data*, JCAP 09 (2024) 016 [arXiv:2404.07203] [INSPIRE].

[39] E. Pajer, D.-G. Wang and B. Zhang, *The UV Sensitivity of Axion Monodromy Inflation*, arXiv:2412.05762 [INSPIRE].
[40] D.-G. Wang and B. Zhang, Bootstrapping the Cosmological Collider with Resonant Features, [arXiv:2505.19066 [INSPIRE]].

[41] N. Arkani-Hamed, D. Baumann, H. Lee and G.L. Pimentel, The Cosmological Bootstrap: Inflationary Correlators from Symmetries and Singularities, JHEP 04 (2020) 105 [arXiv:1811.00024 [INSPIRE]].

[42] G.L. Pimentel and D.-G. Wang, Boostless cosmological collider bootstrap, JHEP 10 (2022) 177 [arXiv:2205.00013 [INSPIRE]].

[43] S. Jazayeri and S. Renaux-Petel, Cosmological bootstrap in slow motion, JHEP 12 (2022) 137 [arXiv:2205.10340 [INSPIRE]].

[44] Z. Qin and Z.-Z. Xianyu, Helical inflation correlators: partial Mellin-Barnes and bootstrap equations, JHEP 04 (2023) 059 [arXiv:2208.13790 [INSPIRE]].

[45] Z. Qin and Z.-Z. Xianyu, Closed-form formulae for inflation correlators, JHEP 07 (2023) 001 [arXiv:2301.07047 [INSPIRE]].

[46] S. Aoki et al., Cosmological correlators with double massive exchanges: bootstrap equation and phenomenology, JHEP 09 (2024) 176 [arXiv:2404.09547 [INSPIRE]].

[47] H. Liu and Z.-Z. Xianyu, Massive Inflationary Amplitudes: Differential Equations and Complete Solutions for General Trees, arXiv:2412.07843 [INSPIRE].

[48] J. Chen, B. Feng and Y.-X. Tao, Multivariate hypergeometric solutions of cosmological (dS) correlators by d log-form differential equations, JHEP 03 (2025) 075 [arXiv:2411.03088 [INSPIRE]].

[49] C. Sleight, A Mellin Space Approach to Cosmological Correlators, JHEP 01 (2020) 090 [arXiv:1906.12302 [INSPIRE]].

[50] C. Sleight and M. Taronna, Bootstrapping Inflationary Correlators in Mellin Space, JHEP 02 (2020) 098 [arXiv:1907.01143 [INSPIRE]].

[51] Z. Qin and Z.-Z. Xianyu, Inflation correlators at the one-loop order: nonanalyticity, factorization, cutting rule, and OPE, JHEP 09 (2023) 116 [arXiv:2304.13295 [INSPIRE]].

[52] Z. Qin and Z.-Z. Xianyu, Nonanalyticity and on-shell factorization of inflation correlators at all loop orders, JHEP 01 (2024) 168 [arXiv:2308.14802 [INSPIRE]].

[53] Z.-Z. Xianyu and J. Zang, Inflation correlators with multiple massive exchanges, JHEP 03 (2024) 070 [arXiv:2309.10849 [INSPIRE]].

[54] Z. Qin, Cosmological correlators at the loop level, JHEP 03 (2025) 051 [arXiv:2411.13636 [INSPIRE]].

[55] Z.-Z. Xianyu and H. Zhang, Bootstrapping one-loop inflation correlators with the spectral decomposition, JHEP 04 (2023) 103 [arXiv:2211.03810 [INSPIRE]].

[56] D. Werth, Spectral representation of cosmological correlators, JHEP 12 (2024) 017 [arXiv:2409.02072 [INSPIRE]].

[57] H. Liu, Z. Qin and Z.-Z. Xianyu, Dispersive bootstrap of massive inflation correlators, JHEP 02 (2025) 101 [arXiv:2407.12299 [INSPIRE]].

[58] J.M. Maldacena and G.L. Pimentel, On graviton non-Gaussianities during inflation, JHEP 09 (2011) 045 [arXiv:1104.2846 [INSPIRE]].

[59] P. Creminelli, Conformal invariance of scalar perturbations in inflation, Phys. Rev. D 85 (2012) 041302 [arXiv:1108.0874 [INSPIRE]].