Analytic Formulae for Inflationary Correlators with Dynamical Mass

Massive fields can imprint unique oscillatory features on primordial correlation functions or inflationary correlators, which is dubbed the cosmological collider signal. In this work, we analytically investigate the effects of a time-dependent mass of a scalar field on inflationary correlators, extending previous numerical studies and implementing techniques developed in the cosmological bootstrap program. The time-dependent mass is in general induced by couplings to the slow-roll inflaton background, with particularly significant effects in the case of non-derivative couplings. By linearly approximating the time dependence, the mode function of the massive scalar is computed analytically, on which we derive analytic formulae for two-, three-, and four-point correlators with the tree-level exchange of the massive scalar. The obtained formulae are utilized to discuss the phenomenological impacts on the power spectrum and bispectrum, and it is found that the scaling behavior of the bispectrum in the squeezed configuration, i.e., the cosmological collider signal, is modified from a time-dependent Boltzmann suppression. By investigating the scaling behavior in detail, we are in principle able to determine the non-derivative couplings between the inflaton and the massive particle.


Introduction
Observations of the Cosmic Microwave Background [1][2][3] strongly support the existence of cosmic inflation [4][5][6][7][8] which can give solutions to several problems in the standard cosmology as well as is the source of primordial scalar (curvature) and tensor perturbations.Furthermore, given the high energy scale (ρ 1/4 inf ∼ 10 15 GeV at most), inflation is a unique opportunity to explore high energy physics including models beyond the Standard Model.
In recent year, there has been a growing interest in an approach that exploit higher-point correlation functions, or non-Gaussianity, of scalar and tensor perturbations, called cosmological collider (CC) program (see earlier works [9][10][11][12] and recent developments ).Indeed, these correlation functions can contain information of (possibly new) massive particles created during inflation, and especially in the soft limits, we expect a sharp oscillatory behavior (what we call "signal") characterized by the mass of particles around the Hubble scale.
A number of recent developments have been seen in the computational methods of primordial correlators.In particular, the so-called cosmological bootstrap method  (see also AdS techniques [139][140][141][142][143][144]) has made it possible to compute the correlation functions rigorously and analytically without having to perform the awkward time integrals of special functions in the cosmological in-in calculations (see Ref. [145,146] for the detail).This allows us to evaluate not only the signal parts of CC but also the non-oscillatory parts (background).This is quite important to understand how large the net signal is.
So far, most of the studies have focused on situations where the masses of massive fields are constant, and consequently, the correlation function has a scale-invariant form.However, interactions with inflaton can produce a non-negligible time dependence on the masses of these massive fields, resulting in a scale-dependent correlation function.For example in Ref. [86], it is numerically shown that a significant deviation from the standard CC signal in bispectrum can be obtained in case of the non-derivative coupling e αϕ/M Pl σ 2 , where ϕ is an inflaton and σ is an isocurvaton.
In this paper, we derive analytic formulae for two-, three-, and four-point inflaton correlation functions or "inflationary correlators".We focus on a simple model consisting of an inflaton and a massive scalar field, where the inflaton imparts a time dependence on the massive scalar through interactions.We then solve the so-called "bootstrap equations" for the correlators [101] by approximating the time-dependence of the mass at the linear order of time.In the constant mass limit, our results consistently reproduce those obtained in Refs.[128,132].The analytic formulae allow us to take into account arbitrary momentum configurations and background parts, so that we can extract various phenomenological aspects of non-Gaussianity in a more precise manner.Specifically, the bispectrum (three-point correlation function of the curvature perturbation) in the squeezed limit contains information on the time dependence of the mass of σ, which is useful to distinguish couplings between inflaton and the massive field.
The paper is organized as follows.In section 2, we introduce the above setup and solve the mode equations for a massive scalar field with time-dependent mass.We also introduce integral formulae for the two-, three-, and four-point correlators of interest.We then derive and solve the bootstrap equations satisfied by the correlators in section 3.In section 4, we discuss the observational impact by taking a particular interaction between the inflaton and the massive scalar.Section 5 is devoted to the summary.Appendices A and B provide boundary conditions necessary to solve the bootstrap equations, and a consistency check in the constant mass limit.

Setup
In this work, we consider a simple system consisting of only an inflaton ϕ and a massive scalar σ, and assume that their interactions include the following non-derivative coupling: During inflation, the interaction (2.1) gives an effective mass of σ, m 2 eff = g(ϕ 0 (t)) with ϕ 0 (t) being the inflaton background, and thus the mass of σ becomes time-dependent.We will first investigate the effects of time dependence with Eq. (2.1) on the mode function of σ.Then, based on the mode function derived there, we study its impact on the inflationary correlation functions in the next sections.

Canonical quantization
We quantize the inflaton fluctuation δϕ and the massive scalar σ with the mode expansion and the commutation relations for the annihilation and the creation operators Note that k expresses three-dimensional vectors, and k ≡ |k| is the absolute value.Under the slow-roll approximation, the mode function of inflaton fluctuations, u k , satisfies an equation of motion for a massless scalar field in de Sitter spacetime, Assuming the Bunch-Davies vacuum, we obtain the canonically normalized mode function On the other hand, the mode function of the massive field, v k , satisfies the following equation: Here m 2 eff ≡ g(ϕ 0 (t)) is an effective mass of σ, which in general depends on time.

Approximation of effective mass for analytic computations
In the forthcoming analysis, we impose several assumptions on the effective mass m eff and employ approximations to perform analytic computations.To explain them, we begin by recalling that, within the framework of the slow-roll approximation, the inflaton background is described by where we used ϕ ′ 0 = √ 2ϵM Pl /τ (assuming φ0 < 0 without loss of generality) and introduced a reference time τ 0 as an integration constant.In this context, the value of the background inflaton field ϕ 0 * at the horizon crossing time τ * for a mode k is given by where we used the condition kτ * = −1 and also introduced a reference scale k 0 such that k 0 τ 0 = −1.The effective mass at the horizon crossing of the mode k reads We are interested in the situation where the variation of the effective mass is not negligible, i.e., |∆m 2 eff | ≳ m 2 eff during observable inflation with e-folding number roughly 50 ≲ N ≲ 60 for CMB, where ∆m 2 eff represents the change in the square of the effective mass from the beginning to the end of observable inflation.It is worth emphasizing that, under this condition, the time dependence of the effective mass produces significant effects that dominate over other slow-roll suppressed effects such as the inflaton mass.Moreover, we mainly focus on the situation that the effective mass remains at the Hubble scale throughout inflation, m eff ∼ H, and undergoes several times larger or smaller changes relative to the initial mass.This assumption is made to avoid exponential suppression of the CC signal ∼ e −O(m eff /H) .
When we solve the equation of motion (2.7) for the massive field σ, it is also necessary to consider the time evolution of the effective mass around the time of the horizon crossing.
In order to estimate the effect, let us perform an expansion of the inflaton background ϕ 0 (τ ) around the horizon crossing time τ * as where the dots stand for higher order terms in 1 − (τ /τ * ).Correspondingly, we expand the effective mass as where g ϕ, * ≡ dg/dϕ| ϕ=ϕ 0 * .The expansion works well at least within a span of several efoldings around the horizon crossing if the following condition is satisfied: In the following analysis, we will take into account the leading order correction, specifically, the second term of Eq. (2.12).On the other hand, the previously mentioned condition |∆m 2 eff | ≳ m 2 eff for non-negligible time dependence during inflation can be rephrased as |g ϕ, * |M Pl /g * ≳ ϵ −1/2 δN −1 with δN ∼ 10 being the change of the e-folding number over observable part of inflation and m 2 eff being evaluated at the horizon crossing time.To sum up, there exists a parameter regime where both conditions are simultaneously satisfied: We work in this regime and study the effects of the time-dependent mass analytically.

Analytical mode function for massive field
With the leading order correction in Eq. (2.12), we simplify the equation of motion (2.7) for the massive field σ to where It is important to note that g * and g ϕ, * are functions of v(k) = k/k 0 , thus µ and κ have scale dependence accordingly.We can solve Eq. (2.15) analytically as where W a,b (•) is the Whittaker function. 2Note that when the mass of σ is constant, g(ϕ) = m 2 0 (and hence κ = 0), Eq. (2.17) is reduced to with µ = (m 0 /H) 2 − 9/4, which reproduces the mode function for a constant mass [9].
• Bulk-to-bulk propagators D ab with a, b = ± for the massive field σ: where θ(•) is a unit step function and

Inflationary correlators
Here, we specify inflationary correlation functions to be investigated in this paper.The correlation functions can be computed from the formula [145], where O is the operator of our interest, and T ( T) is a (anti) time-ordering operator.H I is an interaction Hamiltonian, and all quantities on the right-hand side are understood in the interaction picture.
Figure 1 Diagrams with a massive scalar field σ with time-dependent mass.

Interactions
In addition to Eq. (2.1), we assume the presence of the following derivative (shift-symmetric with respect to δϕ) interactions involving the inflaton fluctuation δϕ and the massive scalar σ, where c 2 and c 3 are coupling constants with mass dimension 1 and −1, respectively.With the presence of the interactions and the formula (2.24), we can calculate inflationary correlators ⟨δϕ • • • δϕ⟩.In particular, this paper focuses on analyzing the two-, three-, and four-point correlation functions of the inflaton fluctuations as shown in Fig. 1.Before proceeding, let us comment on other types of ϕ-σ interactions which arise from the non-derivative coupling (2.1) and have unavoidable contribution to inflaton correlators.The interaction (2.1) can be expanded as and the second term, in conjunction with two transfer vertices (2.25), contributes to the three-point inflaton correlator as a tree, double σ-exchange diagram, in addition to the diagram of our interest (the middle diagram with single σ-exchange of Fig. 1).The contribution from the double-exchange diagram is estimated as while the one from the single-exchange diagram is Here, in order to compare these contributions, let us further assume that the interactions, (2.25) and (2.26), are generated from a single covariant interaction σ(∂ϕ) 2 /Λ with a cutoff scale Λ.Then, the coupling constants are fixed as c 2 ∼ √ ϵM Pl H/Λ and c 3 ∼ 1/Λ, and thus the ratio of these two correlators is evaluated as (2.30) The inequalities are from the lower and the upper bound of Eq. (2.14) with the quasi-single field regime m 2 eff = g(ϕ 0 ) ∼ H 2 .Therefore, we find the contribution from the doubleexchange diagram is comparable to or smaller than the single-exchange one in the cosmological collider setup.In section 4 where some phenomenological aspects of the time-dependent mass are discussed, we choose the model parameters within the region (2.30), in particular near the lower bound.In this case, the contribution from the double-exchange is negligible (about 10% of the single-exchange contribution).Furthermore, each massive propagator is naively expected to be suppressed by the Boltzmann factor e −πµ , leading to the prediction that the double exchange process can receive stronger suppression by one more Boltzmann factor than the single-exchange process3 .Taking the above discussions into account, we ignore the contribution from the double-exchange diagram in this work 4 .

4-point correlators
Let us begin with the correction to the four-point correlation functions since its specific limit produces the three-and two-point correlation functions.
By employing the in-in formula (2.24) (or SK diagrammatic rule [146]), we obtain where the prime on the left-hand side means that a momentum conservation factor (2π is extracted, and "5 per." represents permutations of the external momenta This can be further simplified as where In the first line we inserted the explicit expression of G a , and in the second line we introduced a "seed integral" [128,132], with p 1 and p 2 being constant numbers with p 1,2 > −5/2, and p 12 ≡ p 1 + p 2 .Therefore, once the seed integral (2.33) is computed, the four-point correlation function can be automatically obtained from Eq. (2.32), and this is also true for three-and two-point functions, as shown below.

3-point correlators
In the same way, by utilizing the seed integral (2.33), the three-point function is given by

.34)
The three-point correlator corresponds to a soft limit (k 4 → 0) of the four-point one or the seed integral.

2-point correlators
Finally, the leading order correction to the two-point function is described as follows: Here, we need an expression for double soft limits (k 2 , k 4 → 0) of the seed integral.Note that the above correlator is a correction to the one in free theory,

Bootstrapping Seed Integrals with Time-dependent Mass
Our task is reduced to evaluate the seed integral (2.33), Performing the time-integral in Eq. (3.1) is challenging for general momentum configurations.Therefore, we employ another approach developed in Refs.[101,128,132], where we can obtain analytic expressions for the seed integrals by solving differential equations ("bootstrap equations") for Eq.(3.1).We will see that this method is also applicable to our case.Note that this section describes technical details, so readers who are not interested in the derivation can skip to the final results, Eqs.(3.47), (3.48) for the seed integral, Eqs.(3.53), (3.54) for its single soft limit, and Eqs.(3.59), (3.60) for the double soft limit.

Bootstrap equations
The SK propagators in Eqs.(2.20) and (2.21) satisfy the following differential equations where we recall k s ≡ |k 1 + k 2 |.They satisfy the same differential equations with respect to τ 2 .By introducing with Eqs. (3.2) and (3.3) can be rewritten as where we defined D ab = k 3 s D ab , or explicitly, ) An important observation is that D depends on z i (i = 1, 2) with a specific combination r i z i .Thus, Eqs.(3.6) and (3.7) can be regarded as the differential equations with respect to r i , i.e., and similarly for r 2 .The subsequent task is identifying differential equations (bootstrap equations) for the seed integral.In terms of z i and r i , Eq. (3.1) is expressed as Let us start from the opposite sign seed integral I p 1 p 2 ±∓ .By applying Eq. (3.10), where we used the following formulae for an arbitrary function f (rz), from the third to the fourth line.Thus, we find or equivalently, where Eq. (3.17) gives the bootstrap equations for I p 1 p 2 ±∓ with respect to r 1 .Those with respect to r 2 are obtained in the same way, D p 1 ∓,r 2 I p 1 p 2 ±∓ = 0. Let us shift our focus to the seed integral with the same sign, I p 1 p 2 ±± .The derivation is basically parallel to the previous case, except for the presence of the "source" term on the right-hand side of Eq. (3.7).In fact, we find where the second term of the right-hand side is the new contribution.After performing the one-dimensional integral, we obtain the bootstrap equations (3.20) The same equation for r 2 is provided with a replacement D p 1 ±,r 1 → D p 1 ±,r 2 .As pointed out in Refs.[128,132], it is useful to introduce new variables In terms of u i , the set of bootstrap equations is summarized as and where (3.26)

±∓
The opposite sign seeds I p 1 p 2 ±∓ satisfy the homogeneous differential equations, (3.22) and (3.24).Considering each equation has two independent solutions, a general solution for I p 1 p 2 ±∓ can be expressed as the following combination: where C ±∓|ab are coefficients fixed from boundary conditions later, and U a|b (u) with b = ±1 are the two independent solutions ( Here, we inserted some numerical factors for later convenience, and p F q is defined by where p F q is a generalized hypergeometric function, and the products of gamma function is abbreviated by

±±
The same sign seeds I p 1 p 2 ±± satisfy the inhomogeneous differential equations, (3.23) and (3.25).The general solution is obtained by combining the general solution for the corresponding homogeneous equation with the particular solution for the inhomogeneous one.Therefore, we can take where G p 1 p 2 ±± is a particular solution and the second term is the solution for homogeneous equations consisting of undetermined coefficients C ±±|ab and Eq.(3.28).
Regarding the particular solution, we rewrite the right-hand side of Eq. (3.23) by using where the array on the right-hand side means the binomial coefficient.This motivates us to take the following ansatz for G p 1 p 2 ±± : which can be solved as where (z is the Pochhammer symbol.Therefore, we obtain the particular solution from Eq. (3.37) and Eq.(3.34).While resolving the double summation in Eq. (3.34) is a non-trivial challenge, we manage to perform at least one of the summations, which results in Finally, we note that the particular solution (3.38) automatically satisfies the second inhomogeneous differential equation with respect to u 2 (3.25).This is confirmed by substituting the ansatz (3.34) into Eq.(3.25).Then, in the same way as above, we obtain recurrence relations for χ m,n , The first relation is exactly the same as Eq.(3.35) and one can check the second relation is automatically satisfied by our solution (3.37).Therefore, Eq. (3.38) gives a particular solution for both inhomogeneous differential equations (3.23) and (3.25).

Determining coefficients
We will specify boundary conditions for the seed integrals (3.1) to determine the coefficients C in Eqs.(3.27) and (3.32).In Appendix A, we evaluated the seed integrals in the hierarchical collapsed limit (u 1 ≪ u 2 ≪ 1) by employing another method based on the Mellin-Barnes representation of the Whittaker function [128,132].The results are given by lim where   (3.42).Therefore, the coefficient is determined as C ±∓|ab = C±∓|ab and C ±±|ab = C±±|ab .In the derivation, we used the fact that, in the limit u 1 ≪ u 2 ≪ 1, the particular solution (3.38) is negligible compared to the second term in Eq. (3.32).This is because the particular solution behaves as ∼ u 5+p 12 1 in the collapsed limit, whereas scales for the dominant term is ∼ u

Summary
In summary, the analytic expressions of the seed integrals (3.1) are obtained as follows: where the definition of 2 F 1 is found in Eq. (3.29) and we remind p 12 ≡ p 1 +p 2 and p12 ≡ p 1 −p 2 .These are the generalization from the results for constant mass [128,132] to those with timedependent mass.In fact, they reproduce the results for constant mass when g(ϕ) = m 2 0 , which formally corresponds to take κ = 0 (see Appendix B).
Before going to the three-and two-point correlation functions, let us comment on the momentum dependence of the seed integrals.These are the functions of u 1 and u 2 , and the same is true for constant mass.However, as already mentioned, they also depend on v defined in Eq. (2.9) through µ and κ (see Eq. (2.16)) because of the time-dependent σ-mass.Incorporating this effect into correlation functions is one of the main results of this work, and its impact on cosmological observables will be investigated in the next section.In the following, this additional momentum dependence of the seed integrals is explicitly denoted by I p 1 p 2 ab (u 1 , u 2 , v(k s ))).

Soft limit
The three-point correlation function (2.34) is obtained by taking the soft limit k 4 → 0 (u 2 → 1) of the seed integrals.For the two-point correlation function (2.35), we further take k 2 → 0 (u 1 → 1) corresponding to the double soft limit.These limits are non-trivial at first sight because of the cancellation of some apparent divergences [128,132], so we look at it in detail in the following.

Single soft limit
In the single soft limit k 4 → 0 (u 2 → 1), the opposite sign seed becomes where C±∓|++ is given in Eq. (3.44), and we defined The Fin.{• • • } denotes the finite parts of U p a|b (1), and the divergence is canceled out in the combination U p 2 ∓|+ (1) + U p 2 ∓|− (1).In the same way, for the same sign seed, we have For the first term G p 1 p 2 ±± , only n = 0 of the n-summation in Eq. (3.38) survives in the limit u 2 → 1, and we obtain the expression without the infinite summation, In summary, after some simplification, we obtain in the single soft limit k 4 → 0.Here p12 ≡ p 1 − p 2 and p 12 ≡ p 1 + p 2 .These reproduces the results for constant mass when g(ϕ) = m 2 0 or κ = 0 (see Appendix B).

Double soft limit
To consider the double soft limit, we take a limit k 2 → 0 (u 1 → 1) of the expressions in single soft limit Eqs.(3.49) and (3.51).For the opposite sign seed, we have In the same way as the single soft limit, divergences are canceled in a combination, U p 1 ±|+ (1)+ U p 1 ±|− (1).For the same sign seed, we have where (3.57) Here we used Fin.lim , but with the one from G p 1 p 2 ±± (1, 1, v(k 3 )).In summary, we obtain the double soft limit (k 4 , k 2 → 0) of the seed integrals as As we expected, they reproduce the results for constant mass when g(ϕ) = m 2 0 or κ = 0 (see Appendix B).
To provide a clear and concrete framework for the subsequent discussion, we assume a specific interaction unless explicitly stated otherwise.The first term represents a time-independent bare mass of σ, and the second term introduces an inflaton (or time) dependence.In the expression, α is a time-independent coupling constant, and the limit α → 0 realizes a constant mass scenario.We also note that Eq. (4.1) can be considered as the leading expansion of e αϕ/M Pl , whose impact on the bispectrum was explored numerically in Ref. [86].In this context, the parameters µ and κ in Eq. (2.16) are Furthermore, to avoid a tachyonic mass for σ, we focus on the following parameter region,

Power spectrum
The power spectrum P ζ is defined by where the prime denotes the omission of the momentum conservation factor.The curvature perturbation ζ is related to the inflaton fluctuation δϕ by As is well established, the standard expression for the power spectrum at leading order is given by and scale dependence appears as slow-roll corrections.In our scenario, Eq. (2.35) combined with Eq. (4.5) introduces a correction to the power spectrum, where I p 1 ,p 2 ab with double soft limit are shown in Eqs.(3.59) and (3.60).The correction from a constant mass scalar field was initially computed analytically in Ref. [17] through direct integration.The same outcome was subsequently derived by solving the bootstrap equations in Ref. [128,132].As a consistency check, our result (4.7) reproduces the result in the constant mass limit, i.e., α = 0. Note that, in contrast to the constant mass case where To see these effects, we introduce a ratio, with the denominator representing the correction to the power spectrum in case of a constant mass.In Fig. 2, we plot the variation of R as a function of v(k 1 ).Our choice of parameter includes a fixed m 0 = 2H and ϵ = 0.05/16, with α ranging from −0.6 to 0.6.Depending on the values of v and α, the correction to the power spectrum can either be increased or decreased.Furthermore, we can analytically estimate the leading v-dependence of the power spectrum P (1) where f P (m 0 ) is a function of m 0 .This v dependence coincides with the behavior derived in general single field inflation scenarios (see e.g., [151]).Note that, although the v-dependence is identical, the leading order of slow-roll parameter dependence is different from the general single field case where the leading order is O(ϵ, η).We also note that the v-dependence in Eq. (4.9) is universal as slow-roll correction, while the effect of σ-field on the inflation background is highly model-dependent.For instance, the power spectrum can often be expanded as where ε represents the leading slow-roll order of the theory, e.g., ε = O( √ ϵ) in our case.

Bispectrum
The bispectrum is characterized by the so-called shape function S defined by Combining Eq. (2.34) with Eq. (4.5), we derive the shape function where the single soft limit of I p 1 p 2 ab is described in Eqs.(3.53) and (3.54).In the subsequent subsections, we explore the behavior of the shape function in detail.

Cosmological collider signal at squeezed configuration
Let us consider the configuration k 1 = k 2 , and define x ≡ k 3 /k 1,2 and v ≡ k 1 /k 0 .Consequently, we can express Eq.(4.11) as In this expression, the squeezed limit k 3 ≪ k 1,2 corresponds to x ≪ 1, while the equilateral one The squeezed limit of the bispectrum is of particular interest in cosmological collider physics, as specific oscillatory patterns can sharply appear in this limit.In Fig. 3, we present S/c 2 c 3 √ x as a function of x −1 , with varying values of α.Our parameter choices are v = 1, m 0 = 2H (2.5H) for α > 0(< 0), ϵ = 0.05/16, and M Pl /H = 10 5 . 6It is readily apparent that both the amplitude and the oscillation frequency exhibit strong dependence on k 1 /k 3 .For α > 0, the amplitude experiences damping and the frequency diminish in the squeezed configuration (x ≪ 1).The opposite behavior is observed for α < 0. This feature is more prominent for an increased value of |α|, and the deviation from a constant mass signal is sufficiently significant to observe.Note that these results are consistent with previous works, including calculations within super-horizon approximation [70] and numerical simulation [86].
The physical interpretation of this signal is clear.In the squeezed limit x ≪ 1, we are essentially observing the correlation between a long wavelength mode k 3 and short ones k 1,2 .x as a function of x −1 with several choices of α.Left : Case with α ≥ 0, that is changed as 0 (blue), 0.5 (orange), and 1 (green).Right : Case with α ≤ 0, that is changed as 0 (blue), −0.3 (orange), and −0.5 (green).
The long mode exits the horizon earlier than the short ones and evolves on the super-horizon scale |k 3 τ | ≪ 1, schematically represented by where we showed only τ -dependence of the mode function, disregarding the relative phase of the positive and negative frequency modes.The relative amplitude e −πµ of the positive and negative frequency modes is nothing but (the square root of) the Boltzmann factor responsible for on-shell particle creation relevant to the CC signal.In the present case, µ exhibits scale dependence due to the horizon crossing time difference.For α > 0 (< 0), the mass of σ becomes heavier (lighter) than the constant mass scenario due to its time-dependent nature.Consequently, the mode function is more (less) suppressed by the Boltzmann factor e −πµ compared to the constant mass case, leading to a smaller (larger) amplitude and a shorter (longer) wavelength in the oscillatory signal of the bispectrum.Furthermore, the oscillation frequency is scale-dependent accordingly.

Distinction of couplings with inflaton
Recall that the time-dependent mass appears because the non-derivative couplings break the shift symmetry; it is not easy to obtain similar effects from derivative couplings.Therefore, when we observe the specific behavior in the CC signal, such as the damping/enhancement of amplitude and the shrinking/spreading of wavelength in the squeezed limit, it is attributed to these non-derivative (direct) interactions.Moreover, detailed observations of the timedependent Boltzmann factor associated with non-derivative couplings between the inflaton and the σ-field allow us to distinguish the form of the interactions.This is obtained from distinctive characteristics in the tail of the "scaling" behavior in the CC signal.We will delve deeper into this aspect in the following discussion.
In the shape function (4.12) with Eqs.(3.53) and (3.54), one can extract the Boltzmann factor for µ ≳ 1 with This reveals that the shape function behaves as

S/ √
x ∼ e −πµ(x) x ±iµ + c.c., (4.15) in the squeezed limit, x ≪ 1 7 .Regarding the Boltzmann factor e −πµ(x) , it is crucial to emphasize that µ exhibits x-dependence due to the time-dependent mass arising from the coupling to the inflaton, unlike the constant mass situation 8 .Consequently, the amplitude of the CC signal exhibits a scaling behavior in addition to the well-known oscillation, as observed in the previous subsection (see Fig. 3).In case of the linear coupling (4.1), we can estimate the scaling as a function of x.The left figure of Fig. 4 illustrates the scaling tail of the oscillatory part (orange line) estimated in Eq. (4.16) on the top of the full shape including the background part (4.12).Eq. (4.16) essentially characterizes the scaling behavior in the squeezed limit (x ≪ 1).Remarkably, the scaling behavior exhibits dependence on the coupling between the inflaton and the massive field.For example, in the case of a quadratic coupling, we obtain e −πµ(x) ∼ e −πm 0 /H • e −π m 0 H •βϵ(log vx)(log vx+2) (4.18) instead of Eq. (4.16).This new formulation represents the amplitude of the corresponding CC signal, as demonstrated in the right panel of Fig. 4, and the scaling behavior (x-dependence) differs significantly from the linear case.Therefore, in principle, we can distinguish the coupling g(ϕ) through careful observations and analysis of these scaling behaviors.
It is straightforward to obtain a similar formula for more generic couplings.For example, for a coupling, with n being positive integers, we have the following scaling function, which predicts a different scaling depending on n.

Equilateral configuration
While the CC signal in the squeezed configuration is quite essential to extract information about σ-field, the bispectrum (4.12) exhibits a peak (dominant contribution) at the equilateral limit x ≃ 1.Hence, it is equally valuable to investigate the behavior around the equilateral configuration.In our scenario, the magnitude of the shape function or the non-Gaussianity parameter f NL in the equilateral limit is estimated by using for α ∼ O(1) and m 0 ∼ O(H).
In Fig. 5, we present plots of S for x ≳ 0.1 (equilateral configuration) varying the interaction strength α (left figure) and the additional scale dependence v (right figure).We find that the values of α and v change the scaling behavior of the bispectrum towards the equilateral configuration.Furthermore, the leading v-dependence of the shape function S at x = 1 is obtained in the same way as the power spectrum.Specifically, ∂S/∂v = f S (m 0 ) √ ϵα/v + O(ϵ) where f S (m 0 ) is a function of m 0 .Unfortunately, for reasons discussed in the subsection of the power spectrum, it is difficult to extract information about the theory from the scale dependence at the exact equilateral configuration.

Summary
Our paper provided analytical formulae for the inflationary correlators, specifically two-, three-, and four-point correlation functions, in the presence of a massive scalar field with time-dependent mass.The results are summarized in Eqs.The time dependence of mass naturally arises from couplings to the inflaton, and its effects can be significant when the couplings do not respect the shift symmetry of the inflaton, i.e., in the case of non-derivative couplings.The mode function of the massive field (2.17) is analytically given by the Whittaker function under the linear order expansion around the time of horizon crossing for the time-dependence inflaton background, which characterizes the scale dependence of the signal.The analytical formulae for the correlators obtained from the bootstrap method include both the signal (oscillatory) parts and the background (non-oscillatory) parts of the bispectrum as graphically shown in Fig. 3 and 5.The nonoscillatory part of the signal is also crucial for achieving precision in the cosmological collider (CC) program.
As an application, the phenomenological impact of the time dependence on the primordial curvature power spectrum and bispectrum were investigated.Our analysis revealed that the CC signal exhibits the damping/enhancement of amplitude due to the time-dependent (or scale-dependent) Boltzmann factor.The change in amplitude in the CC signal can be readily observed and provide evidence for the time-dependent mass of the massive field.Although this behavior has been qualitatively confirmed in previous numerical studies such as Ref. [86], we quantitatively parameterized this scaling of the amplitude in a simple form (4.20) based on our analytical formulae.By probing the scaling behavior in the tail of the signal, we are able to distinguish couplings between the inflaton and the massive field, which provides a pathway to explore the unknown inflaton sector through the CC signal.
It should be noted that we made an approximation by linearizing the time dependence of the inflaton and omitting the higher-order time dependence.This is merely the approximation made in our study.While extensions beyond linear order present technical challenges in obtaining analytical representations of non-Gaussianity, there are interesting scenarios beyond the scope of our setup, such as CC signals from tachyonic phase transitions [48].We plan to explore these topics in future work.

±∓
To compute the opposite sign seed integral I p 1 p 2 ±∓ , we first express the propagators D ±∓ using Eq.(A.2): where s 12 = s 1 + s 2 .Substituting this expression into I p 1 p 2 ±∓ enables us to carry out the τ 1,2 -integral, resulting in where and we explicitly denote k s dependence in κ(k s ) (see Eq. (2.16) with Eq. (2.9)).The integration over s i can be computed using the residue theorem.We close the contours from the left with a large semicircle and pick up the poles s i = −n i ± iµ with n i to be non-negative integers.This process allows us to express the seed integral as Let us start from the factorized (F) integral I ±±,F,> .We rewrite the propagators D ±∓ using the Mellin-Barnes representations Eqs.(A.1) and (A.2): which leads to the subsequent expression of I ±±,F,> after the τ i -integration, The integration over s i and the summation of the residues can be executed following a similar procedure, and we obtain where the G-function is defined in Eq. (A.11).
(A.24) Subsequently, the complex s i -integral leads to is obtained by taking a conjugate of Eq. (A.25).In contrast to the previous cases, performing the double summation presented above poses technical challenges.The hierarchical collapsed limit (u 1 ≪ u 2 ≪ 1) of I p 1 p 2 ±± enables us to determine the integration constants appearing in Section 3.Under the limit, the dominant contribution comes from the poles with n 1,2 = 0 in Eqs.(A.21) and (A.25).Additionally, I p 1 p 2 ++,TO,> is subdominant in comparison to I p 1 p 2 ±±,F,> .Combining these evaluations, we obtain lim Note that we confirmed that our outcomes (A.7), (A.22), and (A.25) are consistent with the constant mass results of Refs.[128,132] in the limit κ → 0, as described in the following appendix in detail.

B Constant Mass Limit
In case where the coupling (2.1) is constant, g(ϕ) = m 2 0 , which corresponds to κ → 0, the seed integrals with single soft limit, Eqs.(3.53) and (3.54), should reproduce the results for the constant mass scenario presented in Ref. [128,132].Since the mode function of the massive field σ for the constant mass is expressed by the Hankel function (see Eq. (2.18)) whereas our case is the Whittaker function (see Eq. (2.17)), this limit provides a non-trivial consistency check.
In the subsequent discussion, we frequently use the following formulae for the Gamma function:  These equations are in agreement with the results of Ref. [132].
In the same way, for the seed integrals with double soft limit (Eqs.(3.59) and (3.60)), we can express them in the limit κ → 0 as follows: Again, these expressions are consistent with the results of Ref. [132].

Figure 5
Figure 5 The shape function S as the function of x.Here S is normalized to have value 1 at the equilateral limit x = 1.Left: α-dependence with m 0 = 2H and v = 1 fixed.Right: v-dependence with m 0 = 2H and α = 0.5 fixed.