Cutting Rule for Cosmological Collider Signals: A Bulk Evolution Perspective

We show that the evolution of interacting massive particles in the de Sitter bulk can be understood at leading order as a series of resonant decay and production events. From this perspective, we classify the cosmological collider signals into local and nonlocal categories with drastically different physical origins. This further allows us to derive a cutting rule for efficiently extracting these cosmological collider signals in an analytical fashion. Our cutting rule is a practical way for extracting cosmological collider signals in model building, and can be readily implemented as symbolic computational packages in the future.


Introduction
As the leading paradigm of the primordial universe, inflation provides us more than a natural genesis of a flat, causal, homogeneous and isotropic background spacetime. Not only can the Gaussian quantum fluctuations prepared by inflation seed the large-scale inhomogeneities and anisotropies, but also the deviations from a purely Gaussian spectrum characterize unique physics during the inflationary era. The study of such non-Gaussianities and its implications for the particle spectrum during inflation has become an active field in the recent years, under the title Cosmological Collider (CC) physics . Powered by the fast expansion of the universe, the CC operates at an energy scale as high as H 10 14 GeV, where H is the Hubble parameter during inflation. The inflating spacetime stretches the wavelength of individual modes of a quantum field, leading to the pairwise production of massive particles. If coupled to the inflaton sector, these gravitationally produced particles decay into curvature fluctuations, thereby seeding the non-Gaussian features on the Cosmic that a single tree-level exchange diagram is fully evaluated in a closed form. Yet the bootstrap method highly depends on the dS symmetries, which can be spontaneously broken by a rolling scalar background. For example, in the chemical potential scenario, the rolling inflaton selects a preferred dS patch, enhancing dS-symmetry-breaking particle production. So far, CC signals in this scenario have only been evaluated numerically for particles with mass m 5H, or approximated by using the late-time expansion. However, numerics can be extremely inefficient for large masses, while the late-time expansion does not capture the full waveform of the CC signal, and can even be problematic in certain cases. As a result, how to efficiently extract precise CC signals from the complicated time integrals is a problem faced by all model builders in cosmological collider physics.
We argue that at tree-level, there may be a practical shortcut to solve this problem by considering the physical evolution history in the dS bulk. Although the many nested time integrals obscure the physics, what happens in the bulk at leading order may be described as a combination of three basic events. (i), particles are produced from the background evolution. (ii), particles are produced at interaction vertices. (iii), particles decay into massless inflatons/gravitons. All three types of events happen locally in time, in the form of resonances. Thus the proper time of the particle spent between the resonance events is well defined and corresponds to its dynamical phase. This is the physical origin of the oscillatory CC signals. Hence, depending on the world line of the massive particle, e.g., between (i) and (iii) or between (ii) and (iii), the CC signals can be classified into two categories with different origins and behaviors. On the other hand, these resonances are the stationary points dominating the non-analytic part of the time integral. This motivates us to concentrate on the resonant production/decay picture and propose a practical cutting rule to analytically extract the leading order CC signals for any IR convergent tree diagrams. For the convenience of model builders in cosmological collider physics, we first give a condensed summary of the procedure for this practical cutting rule below and give justifications in the main text.
Cutting rule for CC signals (For a more detailed version, see Sect. 3.3) 1. Given a specific tree diagram, consider only single-colored diagrams with vertices being all-black (+) or all-white (−).
2. Focus on one massive propagator at a time and integrate out all other massive propagators as local EFT operators.
3. Explicitly compute the left/right blobs. The results should be reducible to a sum of polynomials multiplying an exponential function of conformal time.
4. Flip the time-ordering Heaviside step function, and obtain two factorized integrals in the form of Laplace transformations. Discard the leftover commutator integral, which contains no CC signals. 5. Symmetrize and traverse all the left/right injection frequency.

Repeat
Step 2-5 for each massive propagator, then sum them up to obtain the total signal. 7. At last, dress the total signal by the EFT background where all massive fields are integrated out.
This paper is organized as follows. We first classify the oscillatory CC signals and discuss their physical interpretations from a bulk perspective in Sect. 2. Then, in Sect. 3, we introduce the practical cutting rule for the CC signal and justify it using the resonant production/decay picture. This is supplemented with several examples in Sect. 4, including the application to the chemical potential scenario with significant dS symmetry breaking. At last, we conclude and give future prospects in Sect. 5.

Cosmological collider physics in a nutshell
In this section, we begin with a lightning review of cosmological collider physics, and then introduce the resonant production/decay picture for the bulk evolution of massive particles. During inflation, the inflaton background φ 0 usually evolves into a slow-roll attractor phase withφ 0 ∼ const. The scale factor in the FRW ansatz grows exponentially with cosmic time, i.e., a(t) ∝ e Ht = −/Hτ , with H ∼ const being the Hubble parameter and −∞ < τ < 0 being the conformal time. Thus the background spacetime resembles that of dS space in the Poincaré patch. Upon this quasi-dS background, perturbations in the matter sector can be quantized and calculated. In particular, quantum fluctuations of the inflaton field ϕ ≡ φ − φ 0 generate the curvature fluctuation ζ we observed today, via the gauge transformation ζ = −(H/φ 0 )ϕ. These inflaton perturbations may enjoy non-trivial couplings to various massive particles in the dS bulk, and due to their vanishing mass, they survive the dS expansion until the dS boundary at τ = 0. Thus their n-point correlation functions ζ(τ f , x 1 ), . . . ζ(τ f , x n ) | τ f →0 encode interesting information about the particle spectrum of the primordial universe.
To calculate the correlation functions in a given particle physics model with weak couplings, it is most straightforward to use a diagrammatic representation to organize the perturbative expansion. This formalism is known as the Schwinger-Keldysh path integral. The quantum fields start out in the Bunch-Davis (BD) vacuum at τ i → −∞, then evolve according to a local Hamiltonian density, and are finally collected by the n probes placed at τ f → 0. We read out the average value of the n localized probes instead of the probability amplitudes of detecting certain on-shell states which are non-local in spacetime. Due to this special setup, only the "in" state is well-defined but not the "out" state. We also need two copies of evolution histories to accomplish the averaging process, one with time ordering (+) and the other with anti-time ordering (−). More specifically, take the inflaton as the example, the generating functional of non-interacting curvature fluctuations can be written as (2.1) There are totally four types of propagators which can be generated using In momentum space, they are where u k (τ ) is the inflaton mode function. The propagator D ab for massive fields can be derived in a similar way, with u k (τ ) replaced by their own mode functions v k (τ ). As a reference, under the BD initial condition, the scalar mode function reads (2.8) Interaction terms can be perturbatively expanded as usual. After separating the Lagrangian into two parts, the generating functional can be written as (2.10) For weak coupling sizes, we are able to use time-dependent perturbation theory and expand the above functional to any desired order. More specifically, for a given diagram, each internal vertex may be colored to indicate their time ordering. Black vertices (+) represent time ordering while white vertices (−) represent anti-time ordering. For a diagram with V vertices, we must sum the 2 V colored diagrams to obtain the final result. For instance, Fig. 1 shows three typical diagrams one may encounter in cosmological collider physics. As stated in Sect. 1, the reason why correlation functions in dS are difficult to compute is that the time-domain calculation requires doing nested integrals of special functions. Take diagram (b) in Fig. 1, for example, there are effectively 2 6 /2 = 32 colored diagrams for us to calculate . Of these Complex conjugation reduces the number of colored diagram by a factor of two.
-5 -32 diagrams, one contains a 6-layer time integral, four contain a 5-layer integral, and six contain a 4-layer integral, etc. Considering also that the integrand consists of products of ten Hankel/Whittaker functions, it seems rather hard to extract the desired analytical form of CC signals. With the help of effective propagators introduced in [15], one can numerically evaluate this diagram for small masses. However, in the presence of chemical potential, if the mass of the internal particles is large (e.g., m 5H), memory consumption quickly becomes unmanageable due to slow integral convergence. Moreover, the effective propagator is limited to handling two-point mixing between massive field and the external inflaton only, and is not applicable to more general couplings. Therefore, we see that at tree level already, we are in need of a more efficient way to extract CC signals and physics.

Bulk evolution as resonance events
To find an efficient way to extract the CC signals, we should understand the evolution history of massive fields in the dS bulk. First, let us consider a free massive field. Due to the curved background geometry, the mode function of massive field usually takes the form of special functions (e.g. Hankel or Whitaker functions). Yet the mathematical complexity of special functions may conceal the physics inside. Alternatively, a more physical approach is to use the super-adiabatic approximation based on the Stokes-line method [55]. Consider the equation of motion of a canonically normalized massive field with mode function v k I (τ ), The solution in terms of the super-adiabatic approximation reads Here S(τ ) is the Stokes multiplier and F (z) is Dingle's singulant variable, 14) and τ c is the complex turning point in the lower-half complex plane, i.e., w 2 (τ c ) = 0, Im τ c < 0. The above division between the positive frequency part f (τ ) and negative frequency part f * (τ ) minimizes the oscillations in α and β by truncating the super-adiabatic series at an optimal order [59]. This results in the most natural choice of a time-dependent particle number, Therefore, it is clear that gravitational particle production happens when Im F (τ * ) = 0, at Production time: with a duration After production, the initial BD vacuum evolves into a two-particle squeezed state via the Bogoliubov transformation defined by (2.12), For |∆τ * /τ * | 1, which is typical for massive fields, the production of the particle pair can be viewed as an event localized in time. τ * thus marks the starting tick of the clock measured by the dynamical phase of the massive particle. Upon production, the particle pair will propagate in opposite directions and accumulate their dynamical phases in the process. This can be seen more clearly in the Wightman function, (2.21) The late-time limit (τ, τ → 0) momentum dependence defines the local vs non-local separation , The local part is a constant independent of k I and upon Fourier transformation back to space coordinate, gives a contact term proportional to δ 3 (x 1 − x 2 ), hence the locality. The non-local part is non-analytic in k I , with the k 2iµ I term giving rise to long-range correlations of the form |x 1 − x 2 | 2iµ . The non-local part's time dependence manifests the accumulated dynamical phase of two propagating particles, where t stands for the corresponding cosmic time. In contrast, the local part's time dependence is of the form (τ /τ ) ±iµ . This can be viewed as the dynamical phase of one single particle propagating from τ to τ (or τ to τ , but with an extra e −2πµ suppression): Up till now, τ, τ can be any time later than the starting tick at τ * . We will see later that the introduction of interactions automatically leads to events localized in time. These interactions set Another way of this division is in terms of Bessel functions [11], i.e., Jiµ(−kτ ) and J * −iµ (−kτ ). Nevertheless, this division between the local and non-local part is not as sharp. Because before the production time, the negative frequency part has not appeared yet and non-local part propagator should still stay at zero value. the interaction time τ, τ as starting/ending ticks. Before that, we notice a nice property about the non-local part of massive propagator, that it is real after pair production time. (2.25) Consequently, the commutator of field operators does not contain non-local parts: This fact is essentially due to the requirement of microcausality. If the commutator has non-zero non-local parts with k ±iµ I dependence, after Fourier transforming back to position space, the commutator will not vanish out of the light-cone, in direct contradiction with microcausality. The fact above has a direct consequence in the CC signals as we shall see later.
Now let us turn on the interactions of the massive field and the massless inflaton ϕ. Any tree diagram containing a massive field propagator can then be summarized into the form shown in Fig. 2. The left and right blob represent the fully evaluated subdiagram on the left and right, and at time τ, τ , respectively. At early times, τ, τ → −∞, the left (right) blob has a leading time dependence e ik L τ (e ik R τ ), with k L (k R ) containing various momentum combinations from the left (right) subdiagram. For a soft internal momentum k I k L , k R , the massive propagator takes the form (2.22), and the integrals over τ, τ pick up the stationary points These resonances also happen as events localized in time. Physically they can be understood as the decay of the massive particle into massless modes, or its resonant production at the vertex from massless modes. These events effectively mark the ending ticks at τ • , τ • . As a result, we have three types of events happening in the dS bulk. Namely, non-local gravitational production, local vertex resonant production, and local vertex resonant decay. Therefore the massive particle's dynamical phases are accumulated either (I) from the gravitational production event at τ * to the resonant decay event at τ • , τ • , or (II) from the resonance production event at τ • to the resonant decay event at τ • (or vice versa). From this classification criterion, we can deduce two possible types of CC signals, I. Non-local Type:

II. Local Type:
This physical picture of bulk evolution is illustrated in Fig. 2. According to their different origins (dependence about soft momentum), we will name the signals as the non-local type or the local type. Before concluding this section, we make a few remarks on our classification of CC signals.
A more accurate calculation of the stationary points can be obtained via the super-adiabatic mode function, which gives | k 2 Usually, in the soft limit where kI kL, kR, the resonance time reduces to (2.28). Because scale invariant quantities only depend on a subspace in the 3-dimensional parameter space R 3 spanned by k L , k R , k I , namely the projective plane RP 2 . One can freely choose k L k R /k 2 I and k L /k R to parametrize two orthogonal directions in this scale-invariant subspace RP 2 .
2. Physically speaking, these two categories of CC signals have drastically different physical origins. As mentioned above, the non-local signal comes from the gravitational production and vertex decay of two massive particles, whereas the local signal comes from the vertex production and vertex decay of one massive particle.
3. Following the above argument, the two CC signals can have different strengths, because of their intrinsically distinct production mechanisms (i.e., one from background at a linear level, one from perturbation at a non-linear level). The non-local part D non-local −+ is suppressed by |αβ * | ∼ n k I , while the local part D local −+ is, to leading order, independent of the particle number n k I . This is most clearly demonstrated in Sect. 4.3 in the presence of a chemical potential, which can alter the particle number. We will see that the non-local CC signal strength can be greatly influenced by the chemical potential, whereas the local CC signal strength is, to leading order, fixed by dS geometry alone. 4. The resonance picture also has the advantage that it always picks up the leading order contribution, since they represent the stationary points in the time integral. If the stationary point lies in the unphysical region (τ • > 0), there is an overall suppression (−1) iµ = e −πµ , which is the penalty for violating energy conservation in dS.

5.
Notice that in the left and right blobs, the particles need not to be massless. At early times (or equivalently, with hard momenta), all mode functions approach to that of a massless particle with a dynamical phase e −ikτ . Therefore, the intermediate massive particle can resonate with all possible frequencies in the left/right blob, giving rise to various patterns of CC signals . Then the physical interpretation is that the massive particle in consideration can be produced by or decay into all possible particle combinations in the left/right blobs, which are still in the UV "massless" phase.
6. Most of the past studies on CC signals focus on the non-local type, with a few exceptions that explicitly mentioned the local type signal [33,60]. We point out that local signals themselves are not newly revealed result, yet their distinction and interpretation is something we wish to emphasize and exploit in this work.

A Cutting Rule for CC Signals
Based on the above observation on bulk evolution, in this section, we propose a practical cutting rule for the efficient extraction of both types of CC signals. Consider a general n-point tree-level diagram (see Fig. 3) with one massive propagator in the Schwinger-Keldysh formalism. Focusing on the internal massive propagator, the diagram can be expressed as the following form,

The general method
where we have absorbed terms related to the left blob and right blob into F L and F R . In general, they consist of products of the scale factor, coupling constants, as well as the detailed expression of the substructures. We also assume that the time integrals inside the substructures have been performed, so that the only time dependences come from the two endpoints (τ, τ ), to which the In cases with two-point mixing, the early-time limit of the left/right blobs may be in conflict with the late-time limit of the massive propagator. However, one can still obtain the e ikτ dynamical phase by integrating out the massive fields in the left/right blobs as EFT operators of massless fields. In fact, this will be our general strategy to be mentioned in Sect. 3.2. massive propagator is attached. In the early-time limit τ → −∞, or in the soft limit k I → 0 (namely, when the massive mode momentum is much smaller than any momentum combinations in the left/right blob) the time dependence of the left/right blob becomes simple: where k L,R represents various combinations of momentum magnitude sums flowing into the internal massive propagator. Depending on the vertex coloring, the integrand can be separated into the time-ordered part (TO, ab = ++), the anti-time-ordered part (ATO, ab = −−) and the non-time-ordered part (NTO, ab = +− or +−). The NTO part is easier to deal with, because these two time integrals are factorized and can be performed independently. In other words, diagrams with opposite coloring on any adjacent pair of vertices are already cut and become factorized. What actually complicates the calculation is the nested time integral from the TO part , which has bounded integrals such as 0 τ or τ −∞ . However, we find that if one focus on the oscillatory CC signals, which is the most salient feature of cosmological collider physics, then after re-organizing the integrand, one can cut through the time-ordering, and factorize the whole diagram. Consequently, the TO part becomes as simple as the NTO part.
To demonstrate the idea, we focus on the TO part, Without loss of generality, let us assume k L > k R . We can re-organize the integral by flipping one of the Heaviside step function, after which the integration is split into two terms, The first term factorizes into a product of two integrals over (−∞, 0), resembling the easier NTO part. Then we only need to examine the second term with bounded integral, which is essentially the field commutator mentioned in Sect. 2. We will argue that this second term is irrelevant for both types of CC signal if k L > k R . First of all, it is easy to check that this commutator nevertheless cannot contribute to non-local type CC signals. By using the identity (2.26), the commutator from non-local part is exactly canceled for τ > τ > τ * and for τ * > τ > τ . The only non-zero contribution is with the hierarchy τ > τ * > τ . (3.6) The ATO part is obtained from the TO part via complex conjugation and momentum reversal. Recall that the non-local signal requires a pair production event at τ * and two resonant decay events at τ • , τ • . Thus we must require However, (3.6) and (3.7) cannot be satisfied simultaneously, suggesting the absence of the non-local type CC signal. Then the remaining question is whether this commutator contributes to the local type signal. According to our previous physical picture of bulk evolution, the local CC signal can be thought of as the dynamical phase between two resonance events happening at different vertices. The resonance condition at the vertices gives On the other hand, the integration region constrains However, (3.8) and (3.9) cannot be satisfied simultaneously. Therefore, the local type CC signals are also absent in the commutator integral. As summarized in Fig. 4, we have obtained a method to cut the TO integral into factorized integrals, with a small leftover piece free of any desired signals. Note that our initial assumption of k L > k R is crucial for us to argue against the presence of local signals in the leftover commutator integral. However, this is not a limitation to our method. The case with k L < k R is completely analogous. Namely, we only need to flip the other Heaviside step function in (3.4), and the rest of the argument carries over. The complete cut result is then necessarily a piecewise function of k L and k R .
Alternatively, one can obtain the cut result in the parameter region kL > kR and then perform a symmetrization kL ↔ kR.
Before going forward, we can already have a rough estimation of the signal strength using the cut result. After the cut, the signal S > I in the TO part can be estimated as (3.10) Clearly, the two types of signals are suppressed by different factors. The non-local signal is due to gravitational particle production at a linear level, hence it is suppressed by β ∼ n . In contrast, the local signal comes from vertex particle production at a non-linear level, thus it is suppressed by the dS intrinsic Boltzmann factor e −πµ ∼ √ e −m/T dS . This is again the penalty of violating energy conservation in dS as mentioned in Point 4 of the summary above. If the gravitational particle production respects dS symmetries, β = e −πµ and the two signal strengths are degenerate. However, with dS-symmetry breaking chemical potential, the particle number can deviate from that of a simple Boltzmann form, leading to the lifting of the signal strength degeneracy. This is explicitly demonstrated with a vector field later in Sect. 4.3.
As a side remark, it is easy to check the signals in the NTO part, (3.11) Obviously, the non-time-order part is heavily suppressed by |β| 2 , e −πµ |β| or e −2πµ , thus we will henceforth neglect all contributions from NTO parts.

The left/right blob substructures
The above discussion is based on the early-time oscillation behavior of left/right blobs. Namely, the dynamical phase is linearly growing with the conformal time τ . This will indeed be the case if the substructure of left/right blobs contains solely massless fields. Then the massive field will resonate with all possible frequencies contained in F L,R . However, one can imagine that the substructure may also involve massive fields, which allow for dynamical phases growing logarithmically with conformal time. The resonance condition then becomes much more complicated than (3.8), and new types of CC signals would appear. Physically, this corresponds to bulk processes of higher complexity, such as the cascade decay of massive particles. Nevertheless, these processes are heavily suppressed, and there are two ways to treat these extra massive propagators if we focus on the leading CC signal. First, if the massive propagator has nothing to do with the CC signal under consideration, there will be no resonance at the corresponding vertices. Thus it can be integrated out using a partial EFT.
In other words, we mimic its effect by some local operators, whose contribution is perturbative in µ −1 and analytic in the momenta.
Second, if the massive propagator does resonate at the corresponding vertices, the whole diagram will gain more suppression factors (each extra resonant massive propagator contributes an e −πµ suppression). Thus, we can ignore this possibility at leading order safely.
To summarize, our strategy is to work with one resonant massive propagator at a time, while contracting other massive propagators to local EFT operators and absorbing them into the left/right blobs. We then traverse every massive propagator in the tree diagram and sum their CC signals together. This should take care of all leading-order signals non-perturbative in µ −1 . In the end, we can integrate out all massive fields and put back the total EFT tower as a background. Now that we have reduced all massive modes in the left/right blobs, we can work out the possible forms of k L and k R , and hence that of the CC signals. Although their detailed expressions will depend on the substructure of the left/right blobs shown in the Fig.(3), their general forms can be found via mathematical induction. We can prove that in the absence of IR divergences, the time dependence of the left/right blobs is always in the following form: where P L,R are polynomials. The summation on k L,R ranges over all possible frequencies contained in the left/right blob. The proof is quite straightforward. First, notice the integral For a positive integer n ∈ N (which is always the case for IR convergent massless fields), the incomplete gamma function can be simplified recursively to a polynomial multiplied by an exponential function, (3.14) Now, we focus on the substructure of the left blob. We start with the assumption that, a left blob with n L 1 external lines takes the form with P L 1 being a polynomial and k L 1 traversing all frequencies in L 1 . This is obviously true for the simplest blob substructure, namely an n L 1 -point contact vertex. More complex left blobs can be constructed by gluing two smaller left blobs. For instance, consider another left blob L 2 , with Gluing the root of the two trees together by a massless propagator of internal momentum p = whereV 1,2 are interaction vertices that may contain integer multiple of derivative operators and scale factors. The spatial derivatives factor out from the integral, while time derivatives only change the detailed form of the polynomial, which is unimportant here. Hence, we can reduce the blob L 1 ∪ L 2 as follows, (3.18) HereP L 1 ,P L 2 , Q, R, S are polynomials whose detailed forms are of no importance. Notice that in the second step of (3.18), we have formally performed the integral over τ 1 , using the fact that such integrals can be reduced as in (3.14), to polynomials multiplied by an exponential function. Interestingly, after integrating out τ 1 , two different exponential factors appear. The first one depends on the injecting frequency (k L 1 + k L 2 ), the other one depends on the internal momentum (p + k L 2 ). In this way, we have shown that gluing two smaller left blobs gives rise to a larger left blob that shares the same form of time dependence, where the summation ranges over k L 1 ∪L 2 = k L 1 + k L 2 , p + k L 2 for all possible values of k L 1 , k L 2 . By iterating this gluing procedure many times, one can find all possible oscillation frequencies of the left blob L = L 1 ∪ L 2 ∪ · · · . The right blob is treated the same way. Hence all CC signal patterns can be found out as where T denotes all possible combinations of frequency sets.

Summary of the cutting algorithm
In this subsection, we summarize the procedures of applying our cutting rule for CC signals. For a n-point correlation function, denote the set of all massive propagators as Σ.
Cutting algorithm for CC signals 1. Given a specific tree diagram, consider only single-colored diagrams with vertices being all-black (+) or all-white (−). Any mixed-colored diagrams are negligible by at least a factor O(e −2πµ ).
2. Focus on one massive propagator I ∈ Σ and integrate out all other massive propagators as local EFT operators. This reduces the tree topology to that of Fig. 3.
3. Explicitly compute the left/right blobs. The results should be reducible to F L (τ ) = k L P L (τ )e ik L τ and F R (τ ) = k R P R (τ )e ik R τ in the absence of IR divergences, where P L,R are polynomials, and the sums over k L,R take into account of all possible frequencies.
4. For k L > k R , flip θ(τ − τ ) = 1 − θ(τ − τ ) and perform two factorized integrals, The integral takes the form of a Laplace transformation of special functions, which often yields analytical expressions, The leftover commutator integral has no CC signals in it, and can be discarded. 5. Account for the case k L < k R by symmetrization, and traverse all possible frequencies, Before going to explicit examples, we make a few comments on this cutting algorithm.
1. The cut result S + B is only an approximation to the full n-point function. The error comes from two sources. First, to obtain closed expressions, one must truncate the large-mass EFT in both Step 2 and Step 7. This results in errors suppressed by powers of µ −2 , which can be reduced if we include more EFT operators. Second, we have thrown away all the commutator integrals since they do not contain CC signals. However, they are important for restoring differentiability at k L = k R and canceling the folded limit pole at k L,R → k I . This results in an O(e −πµ ) error in these special limits, one that cannot be reduced within our method. Yet in practice, as long as we stay away from these special limits (which only take up a small fraction of the total phase space), this source of error is generically small, as will be demonstrated in Sect. 4.
2. The Laplace transformation of familiar mode functions (Hankel, Whittaker) can be evaluated to hypergeometric functions. This would not have been possible without cutting the nested time integral. In addition, EFT contributions are also analytically calculable. Thus our cutting rule typically results in completely analytical results . Even if the Laplace transform of unfamiliar mode functions may not be evaluated in a closed form, the original two-dimensional integral is now factorized into two one-dimensional integrals. This greatly reduces the computational cost of numerical integration.
3. Each step in the above algorithm is specific but tedious to perform by hand. Therefore, the best way to systematically carry out the cutting procedure is to promote this into a computer program, a goal that we hope to accomplish in the future. This would be beneficial for model builders in cosmological collider physics.
Finally, a comparison must be made between our method and the recently proposed cosmological cutting rules [47,48,[50][51][52]. Based on general principles such as locality, unitarity and analyticity, the cosmological cutting rules give precise relations among the discontinuities of wavefunction exponents Disc k I ψ n . From these relations, one may recursively construct the whole wavefunction. In contrast, our cutting rule directly applies to correlation functions ζ n . It is based on the physical Yet, necessarily non-analytic.
picture of bulk evolution with the specific aim of extracting CC signals of massive fields. Our cut result is also approximate, with reducible EFT truncation errors and irreducible errors from neglecting commutator integrals. However, as stated above, the algorithm itself is well-defined and can be readily implemented as computer programs yielding analytical expressions. We have traded mathematical rigor and formality for practicality and efficiency. Of course, there are also similarities. The fact that the non-local CC signal is associated with a branch cut discontinuity in k I and that it actually comes from a factorized integral seem to hint a connection to the cosmological cutting rules. It is certainly interesting to further explore the similarities and dissimilarities of these two approaches in the future.

Application to Typical Diagrams
In this section, we will explicitly demonstrate how to apply our cutting rule to extract CC signals for the three typical diagrams shown in Fig. 1. For a better illustration, all diagrams calculated in this section correspond to the curvature trispectrum (4-point function). The oft studied bispectrum (3-point function) can be easily obtained by taking the soft limit of one external momentum, and we give more related discussions at the end of Sect. 4.1.

Example 1: Direct scalar exchange
The first example (Fig. 6) is that of a simple exchange diagram with an intermediate massive scalar σ of mass m H = µ 2 + 9 4 . The interaction vertex is chosen to be where ϕ is the massless inflaton. Setp 1, according to the cutting algorithm, we only focus on the TO diagram with the vertices being black. The ATO diagram is obtained via complex conjugation and momenta reversal. More explicitly, Here k I = k 1 + k 2 ,c 3 ≡ c 3φ0 2 H 2 and G, D are the propagators for ϕ, σ, respectively. Without loss of generality, we take the s-channel as an example. The t, u-channels can be included after permutation.
Step 2 and 3 are trivial since no other massive propagators exist in this simple diagram. The left/right blobs are Then at Step 4, denoting k i 1 ···im ≡ k i 1 + · · · + k im , we take k 12 > k 34 and flip one of the Heaviside step functions to split the TO diagram into with the cut result (4.6) and the commutator integral The CC signals can be evaluated via Laplace transformation as Here k 12 is assumed to have a small negative imaginary part in order to select the correct branch of the hypergeometric function. As argued above, the commutator integral contains neither local nor non-local signals. We numerically verify this fact in Appendix A.
Step 5, the total cut CC signals is obtained via symmetrization, Step 6 is trivial and Step 7 is to supplement the CC signals with an EFT background. This background admits a perturbative expansion in powers of µ −1 . These are analytic functions of momenta that can be mimicked by local operators in the single-field EFT. The EFT background The negligible NTO diagram also have a similar form, more details can be found in Appendix B. is usually smaller than CC signals under soft limits. Nevertheless, it plays an important role in "equilateral" regions where all momenta (both external and internal) are comparable in magnitude. In this example, the leading order EFT contribution is generated from the operator In the "s"-channel, its impact on the 4-point function is . (4.11) Finally, we get the analytical approximation for the 4-point function in the s-channel: with S I and B I given by (4.9) and (4.11). In the lower panels, we compare the CC signals filtered from the numerical result and those directly obtained via the cut. We have adopted a high-pass filter with Hamming window.
In Fig. 7, we compare our cut result to that of brute-force numerical integration. The left column shows the non-local type CC signal where we keep k 12 ≈ k 34  In order to extract oscillatory signals from the total numerical result, we have evoked high-pass filters to eliminate the low-frequency EFT background [43]. As we can see from the plots, the numerical results (red dashed lines) match our analytical cut results (blue line) very well, especially for the signal part. The slight mismatch comes around in two places. First, when k 12 ≈ k 34 , the overall signal size differ relatively by O(µ −2 ) ∼ 1 9 , which is the next-order contribution in the EFT tower. We expect that including more EFT operators induced by the exchange diagram can systematically cure this mismatch. The second type of mismatch manifest itself in the folded limit (k 12 , k 34 → k I ) pole as well as in the derivative discontinuity at k 12 = k 34 . As mentioned before, these are caused by the disposal of the commutator integral, which, despite being small, plays an important role in restoring BD vacuum in the early-time limit. This intrinsic mismatch shows a defect of our method. Nevertheless, as long as we stay away from the special configurations, the CC signals we care about will not be influenced much.
Finally, we can gain some insight by studying the large-mass asymptotic behavior of S I in the soft limit, where we have only kept the leading contribution. Clearly, there are two independent oscillation patterns with degenerate amplitudes. In the literature, this large-mass result is often obtained via direct late-time (IR) expansion (2.22). Here, we see that the late-time expansion is qualitatively justified for large masses in this 4-point function case. However, there are two disadvantages for the traditional late-time expansion. First, while the late-time expansion gives the correct exponential suppression factor e −πµ , it may fail to capture the correct mass power dependence (µ n ), which is also important in the estimation of signal strength. Such is the case, for example, in the 3-point diagram frequently encountered in cosmological collider physics. If one performs a late-time expansion before time integration [20], the CC signal amplitude is proportional to µ 2 e −πµ . However, as we can see from tending k 4 → 0 in (4.13), the correct amplitude should be µ 3 e −πµ . This can reproduce the results of [9]. The problem arises because performing a late-time expansion does not commute with completing the time integral. This fact was previously pointed out in [33]. Second, the signal waveform obtained from the late-time expansion is only valid in the soft limit k I → 0. Significant dephasing of the waveform will appear when the internal momentum is less soft. In other words, its validity is limited to a small part of the phase space. In contrast to the late-time expansion, our method can get the correct amplitude as well as the correct waveform of the CC signals, greatly enlarging the phase space we can take advantage of.

Example 2: Scalar exchange with two-point mixing and different masses
The second diagram (Fig. 8) involving two-point mixing is considerably more complicated. As described in Sect. 2.1, the TO diagram alone consists of a 6-layer integral. However, with the cutting rule, we can extract the CC signals without too much effort, as we will demonstrate below. Let the mass of exchanged massive field σ (red line) be µ and the mass of the massive field s on the external lines (cyan) be ν. The interactions are set by the following couplings, (4.14) Without loss of generality, we assume that all propagators are made distinguishable. It is not hard to include permutations according to the Bose statistics. First, focusing on the TO diagram with all vertices being black, the correlation function reads The effective propagator G is defined as [61] G ± (k, where D s and D σ are the massive propagators of the field s and the field σ. Step 2, there are five massive propagators in this diagram. Naming them according to the external lines separated by each one of them, we have the set of all massive propagators, Σ = (12, 34), (1, 234), (2,134), (3,124), (4,123) . (4.17) We focus on one massive propagator I ∈ Σ at a time, integrating out the others as EFT operators. This whole procedure is illustrated in Fig. 9. For example, doing so for I = (12, 34) yields a contact EFT vertex at leading order, ∆L EF T (12,34) = λρ 2 (Hν) 4 a 2 ϕ ϕ σ . (4.18) Step 3, the left/right blobs are now easily computable: (4.20) Step 4, for k 12 > k 34 , we cut the massive propagator I = (12, 34) by with the cut result (4.22) Step 5, symmetrization gives the complete CC signal in the channel I = (12, 34), S (12,34) = θ(k 12 − k 34 )S > (12,34) (k 12 , k 34 ) + θ(k 34 − k 12 )S > (12,34) (k 34 , k 12 ) + c.c. . (4.23) Step 6, we need to loop back to Step 2 and go through the same procedure for the other massive propagators I ∈ Σ. For instance, in the channel I = (1, 234), large-mass EFT expansion of other propagators gives the following local operator at leading order: Thus the left/right blobs are After cut and symmetrization (which is trivial since k 1 < k 234 by the triangle inequality), the corresponding CC signal in the channel I = (12, 34) is then Finally, in Step 7, we dress the signals with an overall EFT background. At leading order, it is sourced by the operator which gives (4.31) As a result, the cutting algorithm gives an analytical approximation to ζ 4 as ζ 4 = S + B, with the detailed expressions of each S I shown in Appendix C. We compare our cut result with numerics in Fig. 10. We have plotted three types of CC signals for two sets of mass parameters. The upper panels correspond to the full 4-point function, while the lower panels show CC signals obtained via the cutting algorithm compared to that filtered from the numeric result. One can see that the oscillatory CC signals match well to that filtered from numerical integration, except at aforementioned special limits where our approximation breaks down. Yet the overall background has an O(1) mismatch. This is because in Step 2 and 7, we have only considered lowest order EFT operators induced by the massive fields . This EFT error is also responsible for the small phase mismatch in the k 234 3k 1 signal. However, we expect such error can be systematically reduced by including higher-order EFT operators. Despite being tedious, the process of including more EFT operators is well-understood and can be performed in an analytical and mechanical fashion. We leave the methodical analysis of the EFT background to future studies.

Example 3: Vector exchange with chemical potential
The third example (Fig. 11) entails a spin-1 massive vector boson with chemical potential. The chemical potential of a vector boson A µ can be introduced via coupling the gradient of a rolling In fact, the EFT error receives amplification from the number of massive propagators. One can estimate the EFT error by R ∼ 2 × 1 µ 2 + 4 ν 2 , which gives R ∼ 30% and R ∼ 52% for the two parameter choices in Fig. 10. This is why it seems to be larger than that in Fig. 7. In the lower panels, we compare the CC signals filtered from the numerical result and that directly obtained via the cut. In cases where the CC signal strength is too weak, we have plotted the analytical cut result only, because the numerical accuracy is inadequate for signal filtering. scalar field θ(τ ) to the Chern-Simons current, The free-theory dynamics of the vector boson is then described by the Lagrangian The mode functions with helicity λ = ±, 0 are expressed in terms of the Whittaker function, where µ = m 2 H 2 − 1 4 ,κ ≡ κ H . The particle production amount can be computed using the Stokes-line method, with the production time and duration given by [55] |k We see that the particle production of transverse modes is suppressed (enhanced) for positive (negative) helicity, whereas the longitudinal mode is not affected by the chemical potential. The vector boson propagators are  After introducing the basic ingredients above, we now compute the 4-point function mediated by the vector boson. We assume the coupling For simplicity, we imagine that the external lines are distinguishable. One can always permute the external momentum after obtaining the result in which they are distinguishable. Following the cutting rule, as a first step, we focus on the TO part of the 4-point function, The kinematic factor reads (4.44) Step 2 being trivial, Step 3 gives the left/right blobs as Step 4, cut the nested integral and discard the leftover commutator integral, where and (4.49) Again, S > I can be analytically computed, we provide its expression in Appendix D . Then in Step 5, we symmetrize the signal and obtain (4.50) Step 6 and 7, the EFT background is sourced at leading order by the local operator In the "s"-channel, its contribution is Finally, we arrive at the result with S I and B I given by (4.50) and (4.52). We compare our cut result to numerics in Fig. 12. Clearly, both the CC signals and the EFT background match the numerical result very well.
To have a closer look at the signal strengths, we take the large-mass limit and the soft limit, then S I becomes (12,34) − e −πµ sin µ ln k 12 k 34 × cos[λφ (12,34) ] . (4.54) We observe that both the amplitude and the phase of the non-local type CC signal is strongly dependent on the chemical potentialκ. For instance, the non-local signal of the negative helicity mode λ = −1 is enhanced forκ > 0, and suppressed forκ < 0. In contrast, its local signal is chemical-potential independent, since it reflects the vertex production/decay process, whose amplitude only depends on the spacetime geometry. This is demonstrated in Fig. 13 by plotting the negative-helicity signals on the projective plane spanned by k 12 k 34 k 2 I and k 12 k 34 . To better visualize the oscillations, we have added a factor (k 12 k 34 /k 2 I ) 9/2 and normalized each plots by their maximum. As we can see from the figure, there are two different oscillations along the horizontal and vertical directions, respectively. In the absence of chemical potential, two CC signals share the same magnitude. However, as soon as we turn on a non-zero chemical potential (|κ| = 0.5), the non-local type becomes dominant (subdominant) if the chemical potential is positive (negative). The positive helicity mode behaves exactly in the opposite fashion, whereas the longitudinal mode is not affected at all. Observationally, the impact of chemical potential on CC signals of different helicities can be distinguished from their dihedral angle dependence. In fact, the case without chemical potential resembles linearly polarized light while that with chemical potential resembles elliptically polarized light. In this way, we can confirm that the two types of signals stem from drastically different phenomena. Notice also that, from the large-mass approximation (4.54), the parity-violating imaginary part of the 4-point function is only present in the non-local signal, not the local signal nor the EFT background. This agrees with the no-go theorem on parity-violating trispectra in single-field EFTs [29]. Figure 13. Signals of the negative helicity component, (k 12 k 34 /k 2 I ) 9/2 S > I (k 12 , k 34 )| λ=−1 , for different chemical potential choices (κ = 0.5, 0, −0.5), with the mass fixed to be µ = 3. Each plot is normalized by their maximum. We see that oscillations in different directions dominates different parameter ranges. Without chemical potential, two CC signals share the same magnitude. As the chemical potential is turned on, the non-local type becomes dominant (subdominant) if the chemical potential is positive (negative).
Before concluding this section, we emphasize that with our cutting algorithm, the numerical efficiency is greatly enhanced. As mentioned in Sect. 1, the conventional numerical integration approach to CC signals is extremely inefficient, and one has to balance between the integration precision and the computational cost. Take the diagram (Fig. 11) considered in this section as an example, we wish to monitor the computational cost of one data point in the configuration space, namely (k 12 = k 34 = 20k I ) for the non-local CC signal of the λ = −1 mode. In Fig. 14, we plot the time (blue line) and memory (red line) consumption for different masses µ, with the chemical potential chosen to beκ = µ − 2. Here we have targeted the relative numerical error at O(10 −3 ). We see clearly that as mass grows larger, both the time and memory cost increase rapidly. For example, when µ = 10,κ = 8, we need around one hour and approximately 3GB of memory to compute even one data point. In contrast, with the help of the cutting rule, evaluating an analytical cut result only costs around O(1)ms and O(10)KB. This is an improvement of 5-6 orders of magnitude, which can be extremely helpful when sampling lots of data points in the trispectrum configuration space.

Summary and outlook
Complicated by nested integrals of special functions, computing inflationary correlators and extracting oscillatory signals are difficult tasks in cosmological collider physics. In most cases, one can only rely on either numerics or work with approximate results of considerable systematic error. In this work, we leverage on the physical picture of bulk evolution, and to make a first step toward a systematic and efficient extraction method of CC signals. We started by a brief review of cosmological collider physics and the Schwinger-Keldysh formalism. After pointing out the technical difficulties for extracting CC signals, we moved on to a physical understanding of the bulk evolution of massive particles. We argued that at leading order, the bulk evolution history can be understood as a collection of resonance events. The oscillatory CC signals can be interpreted as the dynamical phases accumulated between these events, and are classified into two categories with distinct origins. Then, based on this resonance picture, we proposed a cutting rule for extracting analytical CC signals from a general tree diagram. The cutting algorithm is supplemented with three application examples, in which we find good agreement between the cut CC signals and those filtered from numerics. The computational efficiency is also found to be greatly improved compared to traditional numerical methods. As mentioned above, though this work attempts to understand the dynamics of massive fields during inflation and to develop efficient shortcuts to observables, there are certainly many improvements to make and future directions to explore. We conclude by listing some of them below.
• First, our cutting rule is subjected to two sources of errors. The EFT truncation error can sometimes be significant, especially when the particle mass is small or when the number of massive propagators is large. Therefore, we need to take into account the higher order terms in the large-mass EFT in a systematical way. Despite being tedious, this refinement is doable in principle. In contrast, the non-perturbative error at special momentum configurations may require new insights. We hope to reduce both errors in future works.
• Second, the specificity of the cutting algorithm makes it convenient to be implemented as a computer program. We will try to achieve this aim in the future. This may benefit model builders in cosmological collider who wish to skip the technicalities to direct signals.
• Another limitation of our cutting rule is that, so far, it is only applicable to tree diagrams. At loop level, in addition to background pair production and vertex resonant production/decay of massive particles, more complex phenomena may happen, and the cutting rule may need further modifications.
• It is also interesting to study the relation between our phenomenological cutting rule to other more formal cutting rules. As mentioned in Sect. 3.3, they have crucial differences, but their shared features (most importantly, the factorization of non-local CC signals on the branch cut of the internal momentum) seem to hint a hidden connection between them. It is also helpful to understand our physical-picture-based cutting rule from a more fundamental perspective, which may aid to its generalization to other scenarios. 1. The right column shows the absence of local signals in the region k 12 /k 34 > 1. The region k 12 /k 34 < 1 with oscillations is not our concern and is painted gray. Here we have fixed k 12 k 34 /k 2 I = 4 × 10 4 and normalized I I,com by its maximum. 2k I + c.c. .