The cosmological collider in R 2 inflation

Starobinsky's R 2 inflation manifests a best-fit scenario for the power spectrum of primordial density fluctuations. Observables derived from the slow-roll picture of the R 2 model in the Einstein frame relies on the conformal transformation of the metric, which inevitably induces a unique exponential-type couplings of the rolling scalaron with all matter fields during inflation. The “large-field” nature of the R 2 model further invokes non-negligible time and scale dependence to the matter sector through such an exponential coupling, modifying not only the dynamics of matter perturbations on superhorizon scales but also their decay rates. In this work, we identify the simplest observable of the cosmological collider physics built in the background of R 2 inflation, focusing on the so-called “quantum primordial clock” signals created by the non-local propagation of massive scalar perturbations. Our numerical formalism based on the unique conformal coupling can have extended applications to (quasi-)single-field inflationary models with non-trivial couplings to gravity or models that originated from the f(R) modification of gravity.


Introduction
Planck measurements of the primordial scalar power spectrum provide precision tests of the single-field inflationary models within the framework of Einstein gravity [1][2][3].In the latest report [3], the power spectrum of the curvature perturbation with a spectral index n s = 1 has been ruled out by more than 6σ, reinforcing the slow-roll paradigm of singlefield inflation with small deviations from exact scale invariance.Up to 95% confidence level, the observational data favor slow-roll models with concave inflaton potentials, and among representative theories, Starobinsky's R 2 model [4] manifests the best-fit predictions for the central value of the n s distribution and the bounded tensor-to-scalar ratio r.
Next-generation experiments (such as Simons Observatory [5], LiteBIRD [6] and CMB-S4 [7]) will aim to improve the current sensitivity bound on r to be at least one-order-ofmagnitude finer.However, due to the high degeneracy in the predictions from all kinds of viable models, the importance of finding observables to test inflationary theories beyond the joint n s -r constraints has been stressed [8].Recently, a possibility of testing large-field inflation by virtue of the well-established cosmological collider physics [9,50,51] (see also ) has been pointed out in [36].In this work, we further dig into the family of large-field models, focusing on the simplest observable in the cosmological collider for the best-fitting R 2 inflation.This simplest observable, sometimes is referred as the "quantum primordial standard clock [87][88][89]" in the ideal cosmological collider, exhibits the most symbolic prediction from the heavy-mass regime of quasi-single-field inflation [50][51][52][53][54][55][56][57][58][59][60][61][62][63][64].
The ideal cosmological collider is built with approximately exact isometries of the de Sitter spacetime where the analytic structure of correlation functions for the curvature perturbation can be "bootstrapped" from the final boundary surface at the end of inflation without even knowing the detailed time evolution during the slow-roll phase [65][66][67][68][69][70][71][72][73][74][75][76][77][78].In this sense, standard predictions from the cosmological collider physics, typically referring to oscillatory momentum scaling in the spectra of primordial correlation functions created by massive particle productions, are independent of the model that creates the inflationary background.These oscillatory features characterize the non-local propagation of massive particles during inflation, which cannot be mimicked by any effective field theory of the inflaton alone, and they are usually exponentially suppressed and can only exist in three-point or higher-order correlation functions.However, due to the explicit analytic expressions, the standard predictions of the cosmological collider may serve as useful baselines to measure the model-dependent signals from realistic inflationary models.
Realistic slow-roll inflation must at the minimum violate the de Sitter dilational invariance [79] (or the dilatation symmetry) between space and time due to the time evolution of inflaton ϕ = ϕ(t).However, the size of such a dilatation symmetry breaking is suppressed by the smallness of slow-roll parameters and it is in general not easy to result in important observational consequences to the spectra of primordial correlators.For large-field inflation with an inflaton excursion range larger than the Planck scale (∆ϕ ≳ M P ), the pioneer study [36] has identified the effect of dilatation violation in the cosmological collider observables from the scalar particle production with a time-dependent mass.This time-varying effect could be induced by the exponential type coupling ∼ e −αϕ/M P with the slow-rolling inflaton ϕ due to the super-Planckian completion of the inflationary models, and in this case the time-varying effect is only suppressed by the square root of the slow-roll parameter ( √ ϵ ∼ φ/M P H) and it modifies both the spectral amplitude and oscillatory frequencies from the standard predictions with a constant mass.Recently, the analytic expression of primordial spectra featuring such a dynamical mass effect with couplings in the limit of α ≪ 1 is found in [80].
R 2 model is essentially a large-field inflation.However, without additional assumptions, the exponential coupling ∼ e −αϕ/M P with all kinds of matter fields in R 2 inflation naturally arises due to the conformal transformation of the metric into the so-called Einstein frame and α = 2/3 is uniquely fixed by the canonical normalization of inflaton in the Einstein frame.Moreover, such an exponential coupling introduced by the metric conformal transformation (which we shall refer as "conformal coupling" for short) not only acts on the particle masses (or the potentials) but also attaches to the kinetic terms.This can lead to different dilution behavior of the quantum fluctuations in the inflationary background.
The study of the cosmological collider physics in R 2 inflation can have more implications to modified gravity theories beyond the Einstein-Hilbert action and for single-field inflationary models with non-trivial coupling with gravity: • The R 2 model is a special case of the f (R) gravity, while the f (R) gravity can be cast into a generalized version of the Brans-Dicke theory (or more generally the scalar-tensor theory) [81][82][83][84].The Einstein frame of any f (R) gravity theory possesses the same type of conformal coupling with matter fields as studied in this work (see Appendix A.1).  • The Einstein frame of the R 2 model is a useful representation to investigate the slowroll dynamics of inflation.
There has been open debates on the equivalence between the Jordan frame and Einstein frame representations.Since the (conformal) coupling of matter with the slow-roll inflaton (namely the scalaron) does not trivially exist in the Jordan frame, we take the calculation in this work for granted as predictions based on the Einstein frame.
• The inflaton potential of the R 2 model in the Einstein frame coincides with the specific limits of many other inflationary model potentials with more free parameters.Some of those potentials may also be converted via the conformal transformation from their original frame, and this shall generate a similar coupling with matter fields as that in R 2 inflation.In particular, the large-field expansion of the conformal coupling in the Higgs inflation with a non-minimal coupling to gravity [85] takes exactly the same form as that of the R 2 inflation at the leading order (see Appendix A.2).

This work is organized as follows:
In Section 2, we first review the conformal coupling in the Einstein frame of the R 2 inflation and we verify the corresponding interaction vertices induced by such a conformal coupling in Section 2.1.We provide two numerical methods to resolve the mode functions of a scalar perturbation in Section 2.2 & 2.3, and we compare these solutions with an analytic formalism for a time-varying mass in Section 2.4.
We study the leading corrections to the scalaron power spectrum in R 2 inflation generated by the transfer of massive scalar perturbations in Section 3. In Section 4, we clarify the simplest observable of the cosmological collider physics, namely the quantum primordial clock signals from the inflationary background created by the R 2 model.Finally, conclusions and discussions of the main findings in this work are provided in Section 5.

Scalar mode functions
In this work we aim to work out the correlation functions of the curvature perturbation ζ involved with the exchange of a scalar particle σ as illustrated in Figure 1.The physical picture is clear if we restrict ourselves in the spatially flat gauge so that the curvature perturbation is related to the linear inflaton perturbation via ζ = −Hδϕ/ φ0 , where ϕ 0 (t) describes the homogeneous (or the zero-mode) evolution of the inflaton.This simple relation relies on the fact that H and φ0 are nearly constants during inflation.Indeed, this is a good approximation for slow-roll inflationary models with the first slow-roll parameter ϵ 1 ≡ − Ḣ/H 2 ≪ 1. Currently we have a conservative upper bound ϵ 1 < 0.0063 from [65], and the leading corrections to H and φ0 ≈ √ 2ϵ 1 M P H are at least suppressed by factors of O(ϵ 1 ).This means that important time corrections to H or φ0 need to be accumulated over a number of e-folds ∆N > 30 from the pivot scale k * = 0.002 Mpc −1 of observations [65].However, having in mind that the practical experiments could only scan a range of scales across ∆N ∼ O(10), we will treat H and φ0 as constants in time throughout the current investigation.
As we will see, interactions with the inflaton ϕ in the R 2 inflation are typically originated from gravitational effects and thus they all suppressed by the (reduced) Planck scale M P ≈ 2.4 × 10 18 GeV.This ensures that the system considered in this work never enters the "strongcoupling" regime and we can adopt the standard analytic mode functions in the inflationary background for the "free-field" inflaton perturbations.

The R 2 inflationary background
Let us start with a brief review of the slow-roll picture in the R 2 inflation model.More details related to the conformal transformation of the metric can be found in Appendix A. The action relevant to our study is given by: where R J (R E ) is the Ricci scalar in the Jordan (Einstein) frame constructed via the metric ϕ M P makes sure that the inflaton ϕ is canonically normalized in the Einstein frame (2.2).Note that (∂ϕ) 2 E ≡ g µν E ∂ µ ϕ∂ ν ϕ.The mass scale M controls the Hubble scale of inflation H in the Einstein frame and L matter can include all kinds of matter fields during inflation.ϕ is also known as the "scalaron," which can be viewed as the longitudinal degree of freedom of the gravitational sector in the higher-order extension of the Einstein's field equation [81,82].
It is straightforward to compute observables of the R 2 inflation in the Einstein frame based on the scalaron potential (2. 3) The slow-roll parameters of the R 2 model can be obtained from U (ϕ).For example, the first slow-roll parameter reads One can check that the plateau region of the potential U (ϕ) for realizing the slow-roll dynamics (ϵ U ≪ 1) is in the limit of ϕ/M P ≫ 1, which implies that the R 2 model should be cast into the "large-field" inflation scenario.For a given epoch in terms of the e-folding number N ≡ ln a during inflation, the corresponding scalaron value can be obtained via the formula where ∆N = N end − N ≥ 0 is the number of e-folds from the given epoch to the end of inflation.See Figure 8 for a comparsion of (2.5) with the numerical solution of ϕ(N ).This formula allows us to express ϵ U as a function of ∆N .
The most important thing to be outlined here is that: as long as we work out the slow-roll dynamics of the R2 model in the Einstein frame (2.2), the matter sector inevitably receives a universal exponential-type coupling with the inflaton ϕ invoked by the conformal tranformation of the metric.As pointed out in [36], large-field inflation with such an exponential-type coupling can introduce mild, yet non-negligible, scale dependence to the equation of motion of a matter field, modifying solutions of the quantum mode fluctuations from the standard cases with the dilatation symmetry.
To see more explicitly, let us consider as a typical example, L matter = − 1 2 (∂σ) 2 J − V (σ) is a canonical scalar in the Jordan frame (2.1), where (∂σ) 2  J = e ϕ M P g µν E ∂ µ σ∂ ν σ after the conformal transformation.In the Einstein frame (2.2), the equation of motion for the linear perturbation of the isocurvature scalar σ = σ 0 (t) + δσ(t, ⃗ x) becomes where V σσ = ∂ 2 σ V (σ) describes the effective mass around a local minimum at a non-trivial vacuum expectation value (VEV) σ 0 ̸ = 0 such that ∂ σ V | σ=σ 0 = 0. Metric perturbations can be omitted in this equation of motion since the energy density of σ only contributes a negligible fraction to the total.Comparing (2.6) to the case of large-field inflation in [36], one can see that in R 2 inflation the conformal transformation introduces a same time-dependent factor ∼ e −αϕ 0 (t)/M P to the mass term of δσ with α = 2/3.Moreover, the conformal coupling also introduces extra modifications to the Hubble friction term.
If we restrict to the tree-level processes given in Figure 1, the transfer vertices or the interactions vanishes identically from the coupling with potential ∼ e −2αϕ 0 /M P V (σ), since we consider a non-trivial VEV σ 0 near a local minimum with V σ = ∂ σ V = 0.This also makes cubic vertices arising from derivative of the scalar potential, such as V σϕϕ , vanishing.Thus we shall focus on the coupling from the kinetic term.Note that the scalaron ϕ may involve in the scalar potential in the Einstein frame if σ has non-trivial couplings with gravity in the Jordan frame.In such cases σ will become a part of the scalaron under conformal transformation. 2  In the simplest case with a static VEV σ0 = 0, the two-point transfer vertex for the scalar perturbations δσ and δϕ in Figure 1 vanishes from the kinetic term ∼ e −αϕ 0 /M P (∂σ) 2 E and the leading contribution will come from loop diagrams.To allow the presence of the two-point vertex ∼ e −αϕ 0 /M P σ0 δϕδ σ, we assume a zero-mode motion σ0 ̸ = 0 around the pivot scale of our interest, where σ0 must be small enough to preserve a good approximation of the scalaron solution (A.13) in the single-field inflation limit.A possible realization of these conditions is given in Appendix A.3.
In order to examine our numerical computations with the known analytic results in previous studies, we perform integration by parts to the conformal coupling with the kinetic  term in (2.2) to obtain the following vertices and interactions: ) where c 2 , c 3 are constant values of O(1).We will use (2.7) and (2.8) for the calculations of the processes given in Figure 1, which is more convenient for extracting the effects led by the conformal coupling.In the limit of σ0 → 0, the overall tree-level contribution from the original kinetic term ∼ e −αϕ 0 /M P (∂σ) 2 E goes away, which means that the contribution from (2.7) and (2.8) is exactly cancelled out by the other terms from the integration by parts.Thus, with σ0 ̸ = 0, the contribution from the original vertices is expected to be the results based on (2.7) and (2.8) times a suppression factor given by the smallness of the non-vanished σ0 .This suppression factor depends on the explicit model for the isocurvature field.(It is the non-minimal coupling (A.29) if we consider the model of Appendix A.3.) Note that interactions with ϕ originated from the conformal factor are always suppressed by M P , since the scalaron is the longitudinal mode of gravity.This also ensures that mode functions solved by (2.6) are indeed "free-field" solutions in the interaction picture.

The near pivot-scale expansion
We now solve (2.6) by expanding the potential V (σ) around its stable VEV σ = σ 0 to obtain a constant mass parameter V σσ = m 2 σ .Following the approach used in [36], we adopt the expansion around the pivot scale k * for observations: where ϕ * = ϕ(t * ) and a * = a(t * ).As a result, we can rewrite (2.6) into the form of 2ϵ U measures the departure from the exactly dilatation invariant case.Namely, the standard mode function for a massive scalar shall be reproduced in the limit of ∆ ϵ → 0.
The constant ϵ U (and therefore a constant ∆ ϵ ) expansion used in this section allows us to proceed the further transformation: (2.12) Here we consider the canonical quantization of δσ as with In terms of the dimensionless time parameter τ ≡ −kη, which is more convenient for numerical computations, we further translate the equation (2.11) into where M 2 * ≡ m 2 * /H 2 is the dimensionless mass parameter.One can see that the differential equation (2.14) becomes inhomogeneous if ∆ ϵ ̸ = 0, and thus the mode function ũk in general can only be solved numerically.Due to the dilatation symmetry of the de Sitter space, the formulation in terms of the dimensionless time τ = −kη translates the time dependence in the conformal coupling of (2.6) to be explicit k-dependence.
If somehow we could turn off ∆ ϵ only to the mass term in (2.14), the equation of motion becomes This artifical equation exhibits the analytic solution as which recovers the correct Bunch-Davies vacuum state, ũk → −ie iτ / √ 2k, in the early-time limit with τ → ∞.It provides useful information for us to discover the initial conditions for our true solutions. 3et us return to the realistic situation of (2.14).Since we will need to deal with a convolution of time in the correlation functions of the inflaton, for a better numerical performance, we propose a factorization of the mode function as where ν = ν(∆ ϵ , M * ) is defined in (2.16).Such a factorization is motivated by the equationof-motion approach [90] for the strongly coupled regime in the quasi-single-field inflation [59,62,63].Remarkably, we only adopt the factorization to the mode function ũk which has to be solved numerically, while, on the other hand, our system is weakly coupled so that the inflaton mode function takes the standard analytic form of a massless scalar in the de Sitter space.
In terms of the factorized mode function B(k, τ ), we finally arrive at the equation of motion reads A good news is that the mode function ũk (or B) does not recognize the broken dilatation invariance led by a non-zero ∆ ϵ in the flat-space (or namely the early-time) limit with τ → ∞.This means that the mode function shares the standard Bunch-Davies vacuum state, ũk → −ie iτ / √ 2k, with the case of ∆ ϵ = 0.By matching the analytic solution of a standard massive scalar in the vacuum (see Appendix B for more details), we obtain the initial conditions in the limit of τ → ∞ as Here τ UV ≫ 1 represents the UV cutoff in our numerical computation.Some examples of the numerical solutions of the factorized mode function B(τ ) are given in Figure 2. One can see that B essentially captures the oscillation dynamics led by the mass term M * = m * /H in the late-time limit when τ = −kη ≪ 1, while the conventional vacuum mode oscillations (∼ e iτ ) in the early-time limit (τ ≫ 1) has been factored out.The choice of the initial time τ UV ≫ 1 must be sufficiently large such that B can precisely reproduce the standard analytic solution in the scale-invariant case: Note that for the cases of M * > 3/2, one should input ν = iµ with µ = (M 2 * − 9/4) 1/2 to the initial conditions (2.19).In Figure 2, we pick up a constant ϵ U by taking ∆N = 55 in (2.5) as an example so that (

2.21)
There is one more advantage to consider the factorized mode function B. The late-time oscillation of the mode function led by the scalar mass m * is the key to generate the so-call quantum clock signals.For B, this means that the effective mass term in (2.18) shall satisfies or otherwise the mass is tachyonic and the oscillation on superhorizon scales does not occur.
Interestingly, for the presence of a non-zero ∆ ϵ the condition (2.22) becomes time dependent.An illustration of the evolution of the condition (2.22) is given in Figure 3 with the window of constant ∆ ϵ values in R 2 inflation ranging from the choices of ∆N = 50 to ∆N = 60.For the case of ∆ ϵ = ∆ 55 , M * = 1.3 is inside the tachyonic (blue-colored) region around the epoch of horizon crossing at τ = 1.Thus the mode function B is growing on superhorizon scales until sufficiently late time when the scalar mass m * finally dominates and (2.22) is satisfied.The evolution of the M 2 eff (τ ) with M * = 1.3 is shown in the right panel of Figure 3.In this case the mode function has a peak amplitude in the region of τ < 1, as shown in the left panel of Figure 2.This behavior is unfamiliar for the standard case with ∆ ϵ = 0, and it provides a possible realization of the "cosmological tachyon collider [44]."On the other hand, for M 2 * > 2 the condition (2.22) is always satisfied for arbitrary choices of ∆ ϵ .In this case the mode function B starts to decay after crossing the horizon around τ = 1, as shown in the right panel of Figure 2.This is the usual behavior for massive scalars on superhorizon scales.

Solutions with a slow-roll scalaron
For a given potential in slow-roll inflationary models, we can obtain the inflaton field value at the end of inflation, ϕ end , at the epoch when either one of the slow-roll parameters (defined as A.10) reaches to the order of unity.In R 2 inflation, (2.5) gives a good analytic approximation for the inflaton value in terms of the e-folding numbers to the end of inflation (see also where ∆N = ln(a end /a) = ln τ − ln τ end and τ end ≡ −kη end .∆N * = ln τ * − ln τ end is the number of e-folds from k * crosses the horizon to the end of inflation.We take ∆N * = 55 with −k * η * = 1 and −k * η end = e −55 for comparing with results of the pivot-scale expansion approach based on ∆ 55 given by (2.21).Note that τ end ≪ τ IR is required with τ IR being the late-time cutoff in our numerical computations.In terms of the conformal time η, the equation of motion (2.6) becomes where a prime denotes the derivative with respect to η.In terms of the expression in (2.23) for the slow-roll inflaton as a function of the conformal time, we can import the factorized mode function (2.17) and rewrite the equation of motion following the same procedure as where σ and x * = 2/3 ϕ * /M P is the same definition as used in (2.18).The solution depends on, but is not very sensitive to the choice of 50 ≤ ∆N * ≤ 60.Initial conditions of B shares the same as used in (2.19).

Analytic formulae with a dynamical mass
As a reference, we examine the analytic formalism for the scalar mode functions with a timedependent mass provided in [80].In this approach, we simply start with a general assumption for the slow-roll inflation as φ0 = (−Hη)ϕ ′ 0 = √ 2ϵ 1 M P H. Taking ϵ 1 and H as constants, we can work out ϕ 0 as a function of the conformal time η, where ϕ 0 (η) = √ 2ϵ 1 M P ln(η/η 0 ) for some initial epoch η 0 .Expanding the function near a pivot scale k * = −1/η * , we get where ϕ * = ϕ 0 (η * ).Applying this expansion to the conformal coupling with mass, we can obtain where M 2 * = m 2 /H 2 is aligned with solutions in previous methods.For the discussion of R 2 inflation, we impose α = 2/3 and ∆ ϵ = α √ 2ϵ 1 .The equation of motion for the perturbations of a scalar σ with a dynamical mass reads Again, it is convenient to do the transfer δσ = (−kη) 1+∆ϵ/2 δσ and rewrite the equation of motion as where we have defined the parameters

.30)
This equation can be solved analytically with respect to the Bunch-Davies vacuum.The mode function is given by [80]: where W a,b (Z) is the Whittaker W function.
The comparison of numerical solutions from the pivot-scale expansion (Section 2.2) and the slow-roll inflaton approximation (Section 2.3) with analytic mode functions based on the standard constant mass and time-dependent dynamical mass (this section) is given in Figure 4.The standard analytic solution is protected by the dilatation symmetry so that it is invariant under the choice of k.For the choice of k/k * ≫ 1, we find that both numerical solutions agree with each other.

Corrections to the power spectrum
In this section we examine the corrections to the power spectrum P δϕ of ϕ led by the scalar perturbation δσ through the transfer vertex (2.7).This examination can provides a consistency check with the observational constraints [65].In the limit of ∆ ϵ = 0 for the near pivot-scale expansion method (Section 2.2), we can check if the result is also consistent with the closed-form formulae of quasi-single-field inflation presented in [53].
The power spectrum can be obtained from the expectation values of δϕ 2 at the end of inflation through ⟨δϕ ⃗ k 1 )P δϕ , where we adopt the standard analytic form for a massless inflaton mode function as On the other hand, the mode functions of the massive scalar are solved numerically as where ũk = c k Be −ikη and B(k, η) is solved by (2.18) or (2.25).
To compute the correlation functions based on the in-in formalism [91,92], we shall construct the interaction picture of our system.As a first step, we need to identify the Hamiltonian density H[δϕ, δσ, δπ ϕ , δπ σ ] in terms of the conjugate momentum densities δπ ϕ = ∂δL/∂δ φ, δπ σ = ∂δL/∂δ σ, where δL describes the perturbed Lagrangian.At quadratic order of the linear perturbations, the free-field Lagrangian reads This gives the conjugate momentum densities in the interaction picture as After separating the Hamiltonian density into the kinematic part H 0 and the interaction part H I , we replace the perturbations by the ones defined in the interaction picture, namely δϕ I and δσ I , as solutions of the equation of motion given by H 0 .As a final step, we use δ φI = ∂H 0 /∂δπ ϕ and δ σI = ∂H 0 /∂δπ σ to replace the momentum densities so that the Hamiltonian densities are given by We can replace the time-dependent conformal coupling by e − √ 2/3ϕ 0 /M P = e −x * (a/a * ) ∆ϵ for the near pivot-scale expansion approach and by e − √ 2/3ϕ 0 /M P = e −x * (∆N * /∆N ) for the slowroll scalaron approximation, where x * = 2/3 ϕ * /M P .
To compute the n-point correlation function by virtue of the in-in formalism, we perform the expansion: where and the leading corrections are given by the "nontime-ordered" term (3.7) and the "time-ordered" term (3.8).
Taking n = 2 with H I = d 3 ⃗ x H I2 , the correction to the two-point function based on the pivot-scale (PS) expansion approach reads where we have used ϵ U ≈ φ2 0 /(2M 2 P H 2 ) in (3.10) and the numerical factor In the standard case with ∆ ϵ = 0 the factor C PS coincides with C. Taking c 2 = 1 and using ∆ ϵ = ∆ 55 , ϵ U = 2.6 × 10 −4 , ϕ * /M P = 5.26 based on ∆N = 55 in R 2 inflation, we find Here 0|δϕ 2 I |0 ′ = H 2 /2k 3 is the zeroth-order expectation value of the two-point function.
The VEV of σ 0 depends on the model of the isocurvature scalar, yet, in general we expect σ 0 /M P ≪ 1.
Let us focus on the numerical factor C = x 2∆ϵ k C PS in (3.10), which collects the contribution that is independent of the coupling constants of the transfer vertices.In the pivot-scale expansion approach, the factor C(∆ ϵ , M * , x k ) is defined as where such a definition can reproduce the standard analytic results found in [53] when employing the mode function for a massive scalar with iµ (τ ).
(3.13) • Δ ϵ = 0 The numerical factor C(M * ) for the correction to the inflaton power spectrum with k = k * .Results of the near pivot-scale expansion are presented with ∆ ϵ = 0 for the standard dilatation invariant case and ∆ ϵ = ∆ 55 for the R 2 inflation.The scale dependence in the modes functions solved with the slow-roll scalaron cannot be turned off.For a consistency check with the analytic formula (3.17), we provide the numerical results for the standard analytic mode functions using the same set off cutoffs {τ IR , τ UV } for all the numerical computations.
Note that (3.12) is a combined representation for both the time-ordered term (3.8) and the non-time-ordered term (3.7).One can use the trick to the non-time-ordered term and change the order of integration to obtain this represetation.Such a represtation is also used in [51] for the standard ∆ ϵ = 0 case with µ = −iν.It is clear to see that the factor e iτ − e −iτ in the outer integral of (3.12) can cancel out the divergence in the late-time limit with τ → 0.
In terms of the factorized mode function B, the numerical factor is given by for the solution of (2.17).On the other hand, for the slow-roll (SR) scalaron approximation we can follow the same procedure to obtain with the numerical factor defined as The results of C(M * ) computed from different approaches are given in Figure 5.
The standard analytic result based on the dilation invariant mode function (3.13) is found to be of the form [53]: where ψ (1) (z) = d 2 ln Γ(z)/dz 2 is the first derivative of the digamma function.
In Figure 5, we compare numerical results of C(∆ ϵ , M * , x k ) from both methods of the near pivot-scale expansion and the slow-roll scalaron approximation.Their results in R 2 inflation agree with each other.To check the numerical error led by the early-time (or the late-time) cutoff τ UV (τ IR ), we perform the numerical check of (3.14) with B replaced by the standard analytic mode functions (2.20) but using the same set of {τ IR , τ UV }.Since the solution of ∆ ϵ = 0 is supposed to reproduce the standard analytic result, one can see that the error enhances with the increase of M * and therefore we shall truncate the analysis around M * ≲ 6.This numerical error is mainly due to the fast oscillations near the cutoff τ UV , which is independent of our choices of ∆ ϵ or x k .There is an universal enhancement 55 led by the modified decay rate of the scalar perburbation δσ with a non-zero ∆ ϵ .The maximal enhancement is given by where H < 3.04 × 10 13 GeV, or the tensor-to-scalar ratio r < 0.036 is a recently improved upper bound [93].One can check the scale-dependence in terms of is red-tilted and is the same for the two numerical approaches.

The simplest quantum primordial clock
Bispectrum from the 3-point correlation function of the curvature perturbation is the simplest observable for the oscillatory momentum scaling features created by the (virtual) production of a massive particle [9].These oscillatory features are imprints of the quantum interference between vacuum fluctuations of the inflaton (∼ e −ikη ) and the late-time oscillation of the massive particle mode functions on superhorizon scales (∼ e imt ).In the ideal cosmological collider with exact dilatation invariance, such an oscillatory feature is used as an "standard clock" to test the time evolution of the scale factor a(t) in different scenarios for the primordial Universe [87].In this work we only focus on the scenario of inflation, a(t) = e Ht , for the primordial Universe, but we compare the simplest quantum clock signals in R 2 inflation to the prediction in the ideal cosmological collider built with exact dilatation invariance.By "simplest," we mean: (1) It resides in the primordial bispectrum, (2) the generation of such a signal uses the minimal number of interaction vertices, and (3) the shape function of this clock signal can be constructed out of one unified time integration formula. 4s given by the diagrammatic illustration in the right panel of Figure 1, the simplest quantum primordial clock is realized with only two vertices.The two-point transfer vertex given by (3.6) is used in the computation of the power spectrum, and the three-point vertex comes from (2.8) reads Following the in-in formalism, the (simplest) quantum clock signal can be obtained via assuming that δϕ I is a Gaussain random field.
In the computation of this section, we assign k 3 as the external momentum leg fixed by the transfer vertex so that it is conserved with the internal momentum of the massive scalar perturbation where δσ k = δσ k 3 .We always align the equilateral configuration of the kinematic triangle with respect to the pivot scale, such that When the dilatation invariance is broken by the exchange process with the massive scalar perturbation δσ in R 2 inflation, the results from (a) and (b) can be different.We refer the squeezed configuration of (a) as the "forward" scaling since it measures the scales exit the horizon earlier than k * , and thus the configuration of (b) is referred as the "backward" scaling.Note that squeezing the external leg for k 1 or k 2 does not lead to important "clock" signals since the quantum interference only occurs with the mass-induced oscillations in the scalar mode functions u k on scales well outside the horizon.

Shape functions with analytic calibrations
Let us focus on the non-time-ordered term (4.3), which is found to be a subdominant contribution to the "quantum primordial standard clock [87]" in the large mass limit m/H ≫ 1 when comparing to the contribution from the time-ordered part (4.4).However, as we will see, the spectral shape generated by the non-time-ordered term is a pure oscillatory function, Results from the pivot-scale expansion approach S PS are given with the choices of ∆ ϵ = ∆ 55 and ∆ ϵ = 0.The results of ∆ ϵ = 0 in the two panels are in well agreement with the standard analytic prediction and theferore they are identical in the two scaling configurations.We use M * = 2 in this plot.
while the shape function result from the time-ordered integration is in general dominated by some non-oscillatory components.Therefore, the non-time-ordered results can provide a better visualization of the "clock signal."Furthermore, in the standard case with exact symmetries, the non-time-ordered integration can be worked out analytically.This allows us to check if our numerical computation can recover a consistent result by turning off the parameter ∆ ϵ .For the purpose of calibration with the standard analytic formalism, we can only adopt the near pivot-scale expansion approach given by Section 2.2 since the explicit k-dependence in the mass term of the slow-roll scalaron approximation in Section 2.3 cannot be switched off.
The non-time-ordered term (4.3) is composed by the two oppisite time orderings for the interaction vertices given in Figure 1: where these two terms are complex conjugate to each other.To see this, let us first work out the two-point anti-time-ordering and the three-point time-ordering correlation for the pivot-scale expansion approach as being the expectation value of (3.1) at the end of inflation.Thanks to the conformal structure preserved in the inflaton sector (as a massless scalar), it is possible to replace the time evolution of the mode function f k in the bulk by the differential operators defined on the late-time boundary surface according to [9]: where and O B is the boundary operator given by5 Here 2 is used and P , Q are the two independent variables of the kinematic triangle configuration under momentum conservation.They are given by We can now pull the boundary operator out of the time-ordering integral to simplify the calculation.As a result, we find where S PS = x ∆ϵ k 3 S PS = (k 3 /k * ) ∆ϵ S PS collects the spectral shape function that is independent of the coupling constants of the vertices.In terms of the factorized mode function (2.17), the shape function for the simplest quantum primordial clock is constructed by only one integral with its complex conjugation as ) where τ = −k 3 η.For each time-ordering result, S PS is given by Here I PS (∆ ϵ , M * , 1) means the result of the integral (4.12) with P = 1, which corresponds to a "collapsed" triangle of the kinematic configuration.S 32 PS denotes the shape function obtained from the three-point anti-time-ordering and the two-point time-ordering correlation, where one can derive similarly to find that and S 32 PS = x ∆ϵ k 3 S 32 PS .Results of the slow-roll scalaron approximation can be worked out in a same manner.The model-independent shape function is defined according to   The integration of the slow-roll scalaron approximation takes the form of where B is the numerical solution of (2.25).
In the standard case (∆ ϵ = 0) with the analytic solution B(τ )e iτ = √ τ H iµ (τ ) of the scalar mode function, the integrals (4.12) and (4.13) can be resolved analytically.The analytic formula is given in Appendix C. In Figure 6, we present the numerical results of the modelindependent shape function S NT = S 23  NT + S 32 NT from the non-time-ordered term (4.3).The standard analytic predictions (C.1) and (C.2) are in well agreement with the numerical results of ∆ ϵ = 0 from the pivot-scale expansion approach (4.12) and (4. which we can referred as the simplest quantum primordial standard clock.The standard clock signal cares only the shape of the triangle configuration rather than its length k t = k 1 + k 2 + k 3 .However, when the dilatation symmetry is broken by the massive scalar mode function (assigned as k 3 ), the particle exchange process at different scales results in different shape functions.This can be seen from the upper and lower panels in Figure 6 for the shape functions in R 2 inflation scaling over two different regimes.

The clock amplitude
We now include the time-ordered contribution (4.4) to explore the typical size of the entire shape function.The time-order term is involved with two layers of integration.For the pivot-scale expansion appraoch, the time-ordered term is composed as where k 12 = k 1 + k 2 and O B is the boundary operator defined in (4.8).It is convenient to define the time-ordered shape function S TO,PS = x ∆ϵ k 3 S TO,PS with the definition S 32 TO,PS (P ) = where we can reexpress the shape functions S NT,PS = x ∆ϵ k 3 S NT,PS as The structure difference between S NT and S TO given above is very similar to the case of the two-point correlation function.Therefore we perform the same trick as used in (3.12) by decomposing the non-time-ordered integral as After exchanging the order of integration to the second term, we find This result allows us to write down a unified shape function as where S mix is the combined representation of the time-ordered and non-time-ordered contributions given by It is not difficult to see that the combination of e i(P +1)τ − e −i(P −1)τ , with P ≥ 1 being a real number, vanishes in the limit of τ → 0, which regularized the behavior of the shape function in the late-time limit even with ∆ ϵ > 0.
Again, we can repeat the procedure for the slow-roll scalaron approximation to obtain the combined shape functions S mix,SR = x −∆ϵ k S mix,SR : In summary, the three-point correlation function due to the massive scalar exchange process computed in this section reads where S mix = S 23 mix + S 32 mix and S mix,PS = x ∆ϵ k3 S mix,PS for solutions from the pivot-scale expansion approach (2.18) while S mix,SR = x −∆ϵ k3 S mix,SR for solutions from the slow-roll scalaron approximation (2.25).
In Figure 7, we perform a scan of the amplitude for the combined shape function S mix at the equilateral limit P = 2 with k 1 = k 2 = k 3 = k * .The amplitude of the shape function runs from O(0.5) with M * = 2 to O(10 −2 ) with M * = 6.For M * < 6, the numerical errors does not lead to significant deviations to the values obtained from the numerical computations based on the standard analytic mode functions with exact dilatation invariance.
At the end of this section, let us check the typical size of the non-Gaussianity relevant to the simplest quantum primordial clock.We adopt the definition of the shape function from the curvature perturbation as [92]: where (2π is the amplitude of the power spectrum of the curvature perturbation ζ.This gives the definition of f NL in the equilateral limit as where we have used , and the result (4.37) for both numerical solutions lead to the similar estimation where the maximal value (k * /H) −∆ 55 ≈ 10 is given by the highest inflationary Hubble scale H = 3 × 10 13 GeV according to (3.18).The amplitude of S eq mix is shown in Figure 7.

Conclusions and discussions
In this work we have investigated numerical solutions of the quantum mode functions of a massive scalar field during R 2 inflation, focusing on the time-varying mass effect led by the conformal coupling with the rolling scalaron.We have examined the tree-level corrections to the primordial power spectrum and the simplest realization of the so-called quantum primordial clock signals induced by the non-local propagation of the massive scalar perturbations under the slow-roll background created by the R 2 model.These tree-level corrections can occur with a small but non-vanishing homogeneous (or the zero-mode) motion of the massive scalar near a local potential minimum.As shown in Appendix A.3, the isocurvature scalar motion can be supported by a non-minimal coupling with gravity while keeping a negligible backreaction to the scalaron motion along the standard R 2 inflationary trajectory.The most important findings of this work can be outlined as: • The slow-roll inflaton ϕ (or namely the scalaron) in the Einstein frame of R 2 inflation couples to all matter fields through the unique factor e −αϕ/M P with α = 2/3 under the conformal transformation, leading to the breaking of dilatation symmetry in the mode functions of matter perturbations.For a massive scalar field σ, the slow-roll scalaron results in a time-varying mass m 2 (t) ∼ e −αϕ(t)/M P m 2 σ and a modified decay rate δσ ∼ a −(3+∆ϵ)/2 on superhorizon scales, where ∆ ϵ ≈ 0.018 (2.21).
• For a pivot scale k * = 0.002 Mpc −1 , the modified decay rate for a massive scalar perturbation δσ results in a universal enhancement (k * /H) −∆ϵ ≲ 10 to all primordial correlators involved with the exchange process for δσ.The maximal enhancement is given by the highest inflationary Hubble scale H = 3 × 10 13 GeV (3.18) from the recent bound [93].For the corrections to the primordial power spectrum, we find where β model ∼ 10 −8 × σ 2 0 /M 2 P depends on the models of the massive scalar σ, σ 0 is the scalar VEV, and P 0 is the scalaron power spectrum without corrections.∆P is the contribution from the vertex (2.7) only.The numerical factor C is given by (3.14) or (3.16).
• The explicit scale and time dependence in the mode functions of the scalar perturbations brings in scale dependence to the shape functions of primordial bispectra.Figure 6 shows the discrepancy of the spectral shape functions in the same (equilateral-to-squeezed) configuration but scaling over two different ranges of k-scales.The typical size of non-Gaussianity corresponding to the exchange process led by (2.7) and (2.8) is with the same β model from the power spectrum and 2ℜ[S eq mix ] ≲ 0.5 is shown in Figure 7.Given that we only consider vertices arisen from the conformal factor via integration by parts, each vertex of (2.7) or (2.8) is at least suppressed by M −2 P .This already makes all the relevant predictions very difficult to be tested by future experiments.Futhermore, the total tree-level contribution in the standard R 2 inflationary background is suppressed by the smallness of the model-dependent zero-mode motion of the isocurvature scalar, and in general we expect leading corrections come from loops if the scalar is completely frozen at a local minimum.Nevertheless, the investigation is still worth an effort, since R 2 inflation is the best-fit scenario.Moreover, we have not yet explored the full possibilities to obtain some much larger β model from different model buildings.Following the findings of this work, there are several topics remaining to be studied.For example: • Higgs inflation.In the large-field regime (ϕ/M P ≫ 1), the conformal factor in the Higgs inflation with a non-minimal coupling to gravity [85] coincides with that in the R 2 model (see Appendix A.2).In Higgs inflation, one can write down direct couplings of Higgs with scalar fields, for example, in the form of h 2 σ 2 ∼ M 2 P ξ e αϕ/M P σ 2 .Direct couplings of the inflaton with scalar fields can lead to unsuppressed interaction vertices.Of course, the time dependence of the scalar mass could be changed drastically by the introduction of such a direct coupling.
• Particles with spins.The characteristic signals for spin 1/2 fermions and gauge bosons are also important targets for the cosmological collider program [27,28].The Lagrangian which describes the dynamics of these matter fields with non-zero spins must include covariant derivatives with respect to the spacetime metric.Moving to the Einstein frame for the inflaton via the conformal transformation, those terms involved with covariant derivatives will partially cancel out the conformal factor Ω, leading to a breaking of the dilatation invariance in the equation of motion as the case for a real scalar field investigated in this work.
• The Standard Model mass spectrum.The conformal coupling acts on all matter fields, including the Standard Model particles.The clarification of the mass spectrum of the Standard Model is the important first step towards the discovery of new particles by using the cosmological collider [15,16].If R 2 inflation is the correct model, the Standard Model predictions shall be modified by the effects considered in this work.
We leave the further researches of these interesting topics as future efforts.
with Ω(ϕ) = e ϕ/( √ 6M P ) .In terms of the dimensionless parameter x ≡ 2 In Figure 8, we compare this analytic approximation with the numerical solution of the scalaron dynamics ϕ(∆N ) in the potential (2.3) without imposing any slow-roll conditions.Assuming that the slow-roll phase ends around ϵ U = 1, which indicates ϕ end /M P ∼ 1, we confirm that (A.13) is a good approximation for the epoch of inflation (−60 < N < −40) considered in this work.

A.2 Higgs inflation with a non-minimal coupling to gravity
Let us also take a look at the inflationary model driven by the Standard Model Higgs boson with a non-minimal coupling ξ ≥ 0 to the Ricci scalar [85].In the Jordan frame, the action of this model is given by where h is the real scalar mode of the Higgs field in the unitary gauge and we drop the electroweak vacuum expectation value since it is irrelevant to our discussion.Similarly, we shall assign the conformal factor with respect to the scalar function to obtain the Einstein frame representation of the theory (A.14), which is given by Note that the total derivative term in (A.2) has been dropped.The canonical inflaton in the Einstein frame can be derived from [85]: . (A.17)

Figure 1 .
Figure 1.Diagrammatic illustration of the correction to the power spectrum (Left Panel) computed in Section 3 and the corresponding three-point correlation for the simplest quantum primordial clock (Right Panel) in Section 4.
E µν and the definition of the conformal factor Ω 2 = e 2 3
)where t * (or η * ) is the (conformal) time when the pivot scale k * crosses the horizon at which −k * η * = 1.The valid range for this expansion is k min ≪ k ≪ k max , with k min = e −10 k * and k max = e 10 k * .Note that φ0 < 0 in R 2 inflation.As we get rid of the mild time dependence in H and φ0 , we obtain a constant slow-roll parameter ϵ U (t * ) = φ2 0 (t * )/2M 2 P H 2 (t * ).Now the non-negligible time dependence shows only in the exponential conformal factor, as e

Figure 4 .
Figure 4. Numerical solutions of the factorized mode function B with k/k * = 1 (Left Panel) and k/k * = 1000 (Right Panel) based on the pivot-scale expansion (Section 2.2) and the slow-roll scalaron approximation (Section 2.3) with M * = 2 and ∆ ϵ = ∆ 55 given by(2.21).These results are compared with analytic formulae with a standard constant mass and a dynamical time-dependent mass (Section 2.4).

Figure 6 .
Figure 6.The model-independent shape function from the non-time-ordered contribution (4.3) for (Upper Panel) the forward-scaling triangle configuration with k 1 = k 2 = k * , k 3 ≤ k * and (Lower Panel) the backward-scaling triangle configuration withk 1 = k 2 ≥ k * , k 3 = k * .Results from the pivot-scale expansion approach S PS are given with the choices of ∆ ϵ = ∆ 55 and ∆ ϵ = 0.The results of ∆ ϵ = 0 in the two panels are in well agreement with the standard analytic prediction and theferore they are identical in the two scaling configurations.We use M * = 2 in this plot.

3 ϕMe
P , one can obtain the slow-roll parameters from the potential −x (2e−x − 1) (1 − e −x ) 2 .(A.11)For a given value of ϕ, the corresponding e-folding number N = ln a to the end of inflation can be estimated from[94]:∆N ≡ |N − N end | ≈ 2B * (τ 2 ), (4.26) and B is the numerical solution of (2.18).As a result, we can rewrite (4.24) as