Precision Calculation of Inflation Correlators at One Loop

We initiate a systematic study of precision calculation of the inflation correlators at the 1-loop level, starting in this paper with bosonic 1-loop bispectrum with chemical-potential enhancement. Such 1-loop processes could lead to important cosmological collider observables but are notoriously difficult to compute due to the lack of symmetries. We attack the problem from a direct numerical approach based on the real-time Schwinger-Keldysh formalism and show full numerical results for arbitrary kinematics containing both the oscillatory"signals"and the"backgrounds". Our results show that, while the non-oscillatory part can be one to two orders of magnitude larger, the oscillatory signal can be separated out by applying appropriate high-pass filters. We have also compared the result with analytic estimates typically adopted in the literature. While the amplitude is comparable, there is a non-negligible deviation in the frequency of the oscillatory part away from the extreme squeezed limit.


Introduction
Cosmic inflation is currently the leading paradigm explaining the origin of the large-scale inhomogeneity and anisotropy of our universe.The inflation is believed to take place at a very high energy scale, up to 10 14 GeV in terms of the Hubble parameter H.At such high energies, quantum fields, including the spacetime itself, experience strong quantum fluctuations, and these fluctuations can imprint the large-scale inhomogeneity by coupling to the spacetime curvature perturbation.The n-point correlation of the curvature perturbation, or the primordial non-Gaussianity as it is often called, can then record the very high energy dynamics that happen during inflation [1].
In particular, it has been suggested recently that the soft limits of the n-point correlators can be the discovery channels for heavy particles and new interactions with masses up to O(H).These studies are based on earlier works on primordial non-Gaussianities and are dubbed "cosmological collider physics."[2][3][4][5][6][7] The general idea is that a heavy particle with m ∼ H can be created from the vacuum quantum fluctuation during inflation and its physical momentum then quickly redshifts to essentially zero.Being a non-relativistic state, its wave function would oscillate with a fixed physical frequency m, and this oscillation can interfere with the mode function of the curvature perturbation ζ, producing a characteristic oscillatory signal in various soft limits of ζcorrelators, including the "squeezed limit" (k 3 k 1 k 2 where k i = |k i |) of the 3-point function (bispectrum), and the "collapsed limit" (|k ) of the 4-point function (trispectrum).Later on, this program is generalized to include new mechanisms of generating primordial fluctuations and new ways of producing heavy particles during inflation, alleviating several constraints in the vanilla slow-roll inflations and also producing particles with masses much higher than H [8][9][10][11].At present we have already quite a few particle physics models that are capable of generating visibly large cosmological collider (CC) signals .These signals could be searched for in the large-scale structure surveys in the near future or more futuristic 21 cm tomography from the dark ages [33][34][35].
Considerations from particle physics model building show that the 1-loop process could be important for CC signals.For example, a promising class of signals comes from the chemicalpotential-enhanced particles with nonzero spins, including both fermions and gauge bosons [9,10].For this class of models, the chemical potential is provided by the rolling inflaton via the axiontype dim-5 coupling to the fermions or gauge bosons [8,[36][37][38][39][40].Such enhanced states are always in transverse polarizations and thus can enter the 3-point function only through loops.Analytical estimates showed that such loop correlators can be quite large and give rise to large CC signals [9,10].Therefore at least in this type of models, the loop process is the leading contribution to the signal while the tree-level processes are absent.It is therefore important to have a controlled and reliable procedure for computing these loop correlators.
It turns out that the computation of cosmic correlators is highly nontrivial and challenging.There is a well-developed in-in formalism that essentially reduces the cosmic correlators to Feynman diagrams, but with several complications compared with their flat-space counterparts [18,41].First, due to a lack of symmetry, space and time cannot be dealt with at the same footing, which renders the usual 4-momentum representation unusable.Second, the computation of correlation functions involves complicated time ordering, the complication of which grows fast with the number of interaction vertices.Third, the mode functions of particles in inflation involve products of special functions that are rather intractable analytically.Essentially due to these complications, the computation of cosmic correlators is still in a primitive stage.Compared to the highly developed and almost industrialized Feynman diagram computation in flat space up to rather high loop orders, a systematic computation of cosmic correlators of simplest processes at tree level are only achieved recently, while the full 1-loop results are known only for a handful of simple diagrams.
The recent progress [42][43][44] in the analytical computation of tree-level diagrams is essentially achieved by fully exploiting the symmetry of the inflation background, which can sometimes be approximated by de Sitter (dS) space.Due to the enlarged symmetry of dS relative to a general slow-roll background, it is possible to transform a flat-space correlator into the corresponding object in dS via dS time translation and boosts, a process called "cosmological bootstrap".This technique is very useful for the process that respects the full dS symmetry.An alternative approach based on AdS techniques can also handle these dS invariant processes [45,46].There is another line of research aiming at understanding the analytical structure of cosmic correlators in recent years, or trying to go beyond the full dS covariance [47][48][49][50][51][52][53].
Compared with the tree-level progress mentioned above, much less is known about the 1-loop processes.Even worse, the phenomenologically important loop processes with chemical potential enhancement do not respect the full dS symmetry, which makes the aforementioned symmetrybased techniques not directly applicable.Our current understanding of these loop processes is thus only based on analytical approximations that might be applicable in special limits, but are not fully justified in more general cases.In view of this situation, it is desirable to have a fully controlled numerical computation for such processes.
In this paper, we attempt a full numerical implementation of such 1-loop diagrams at the 3-point level.We use the standard Schwinger-Keldysh (SK) formalism in real-time, which does not rely on the dS symmetries.We use the diagrammatic representation and follow the Feynman rules presented in [18], and then numerically carry out all the integrals.In particular, we do not expand the mode function around any particular point as was often done in previous analytical estimates.Instead, we adopt a piecewise expansion of the mode function, making sure that the mode function being used agrees with the full result up to controllable numerical errors.
As far as we know this is the first full numerical computation of 1-loop cosmic correlators that are relevant to cosmological collider physics.The signals from these processes have been studied analytically with several approximations made, some of which are valid only in the squeezed limit of the bispectrum and are applicable only to the oscillatory part of the process [9,10].The numerical approach we adopt here is mostly free of those approximations and thus is valid for arbitrary configurations.The numerical result we obtain can also cover both the "background" and the "signal" part of the bispectrum, and thus is useful for generating templates when confronting the theoretical predictions with data.The "signal" part obtained in this approach is also a useful check of the analytical estimates used in previous works.
As the first step of this numerical program, we consider bosonic processes in this paper, including the loop diagrams mediated by massive scalars and gauge bosons.We include in particular the axion-type dim-5 couplings φF F to introduce a chemical potential enhancement.There are subtleties peculiar to the numerical calculation which we will explain in full detail.Our results show that the overall amplitude can be one or two orders of magnitude larger than the oscillatory signal.At the same time, when the oscillation frequency is not too small, the oscillatory part of the signal can be filtered out from the full result by simple filtering techniques.We show that the oscillatory signals have overall amplitudes and scaling behavior consistent with analytical estimates.At the same time, there is a non-negligible deviation in the frequency of the oscillatory part for k 1 /k 3 < 20.
The rest of the paper is organized as follows.In Sec. 2 we briefly review the cosmological collider physics and its observables.We first introduce the signal from a purely phenomenological point of view and then introduce the diagrammatic formalism for the actual computation.In Sec. 3 we introduce particle physics models that can give rise to CC signals at 1-loop order, including the case of a massive scalar loop and a massive gauge boson loop, both with or without chemical potential enhancement.Then in Sec. 4 we explain in detail the numerical calculation of these 1loop correlators and present our results.Further discussions are collected in Sec. 5. We crosscheck some key aspects of numerical implementations using tree-level correlators in App.A, discuss the behavior of the loop integrand at the large loop momentum limit in App.B, and crosscheck the filtering methods we adopt to separate the oscillatory and non-oscillatory parts of the signal in App. C.

Brief Review of Cosmological Collider Observables
In this section, we review the basics of the CC observables.The program of cosmological collider physics focuses on the inflaton-spectator interaction, aiming to extract the information of those spectator fields from the inflaton correlators.This is interesting in particular because the spectator fields could have too large a mass to be reached by any terrestrial experiments, but are not much heavier than the Hubble scale of the inflation so that they can be effectively produced during inflation.

Observables
The CC signals appear in various soft limits of n-point cosmic correlators.Among them, the squeezed limit of the 3-point correlator is the simplest and therefore will be our main focus.By the translation and rotation symmetries of 3-space, the 3-point correlator is a function of the triangle formed by the three external 3-momenta k i (i = 1, 2, 3) and it depends on k i only through their magnitudes, k i ≡ |k i |.The squeezed limit then refers to the limit when one of k i is much smaller than the other two, e.g., k 3 k 1 k 2 .For the 3-point correlator, it is conventional to define a dimensionless shape function S(k 1 , k 2 , k 3 ) in the following way, where a prime in • • • means to strip away the momentum-conserving factor (2π) 3 δ (3) (k 1 +k 2 +k 3 ), and P ζ 2×10 −9 is the nearly scale-invariant power spectrum, which measures (the square of) the size of the curvature perturbation ζ.In the single-field slow-roll models, ζ is related to the inflaton fluctuation ϕ via ζ k −(H/ φ0 )ϕ k , where φ0 is the rolling speed of the inflaton background.The approximate scale invariance dictates that the shape function S(k 1 , k 2 , k 3 ) depends only on k-ratios up to slow-roll corrections so that it is really a function of "shape" rather than the "size" of the momentum triangle.Consequently, there are only two independent variables that the shape function can depend on.In the squeezed limit k 3 /k 1 → 0 and k 1 k 2 , it is convenient to choose the ratio = k 3 /k 1 and the angle ϑ between k 1 and k 3 as the two independent variables, so that S = S( , ϑ).The ϑ-dependence contains the information about the angular momentum of the particles mediating the process, including both the intrinsic spin and the extrinsic one, and this angular dependence is completely fixed by angular momentum conservation.The dependence, on the other hand, contains interesting information about the mass and possibly other dynamical properties of the intermediate particles.
In fact, in the presence of a heavy intermediate particle, the shape function is generally not analytic in the → 0 limit but can contain a nonanalytic piece in the form of a noninteger power α where α is generally a complex number.Therefore, we expect that the shape function behaves like (suppressing the angular dependence) Here we have several parameters, among which A and B are real numbers measuring the sizes of analytic and nonanalytic pieces, respectively.In general, A can be either positive or negative, while we can always choose B to be positive.The exponent N is an integer, while L is in general a real number.When ω = 0, the exponent L is often an integer or half-integer, although exceptions exist [54].In the nonanalytic piece, we see the characteristic oscillatory dependence on , described by two parameters, the frequency ω and the phase ϕ.
It has been noticed that the dimensionless frequency ω can often be related to the mass of the intermediate particle.For example, the tree-level exchange of a scalar particle of mass m gives ω = m 2 /H 2 − 9/4 when m > 3H/2.But here we emphasize that all four parameters (B, L, ω, ϕ) can in principle be unambiguously measured, and they can provide very useful information about the dynamics of the process.Therefore it is important to calculate all these parameters from a given model.Below we discuss the analytic and nonanalytic pieces respectively.Analytic Piece.In general, the analytic piece A N can be thought of as a result of local (pointlike) interactions.The local interaction can be either from the self-interaction of the external modes or from "integrating out" the intermediate particle when it is heavy (m H).So it corresponds to the effective field theory (EFT) limit when m H. Therefore, in the following we shall also call it the "EFT piece," although this piece is still present even when the intermediate particle is not so heavy (m ∼ or < H) and the EFT limit does not make clear sense.
The integer power N in the EFT piece contains information about the EFT coupling of the external modes.For instance, N = 1 when the external modes are derivatively coupled and N = −1 when they are gravitationally coupled.Direct (non-derivative) coupling of external modes also gives N = −1.But in this case, the amplitude A would have a mild (logarithmic) dependence on since the direct coupling softly breaks the scale invariance in the infrared (IR).
Let us also comment on the terminology which could be rather confusing to nonexperts.We already mentioned that the shape function is a function of the shape of triangles.In literature, a particular dependence on the shape is often simply called a "shape."Often appeared shapes include the local shape and the equilateral shape.Both of them are relevant to our study of cosmological collider physics and we will comment on them below.Rather confusingly, people also occasionally talked about the shape of the momentum triangle, such as the equilateral shape k 1 = k 2 = k 3 , the squeezed shape k 1 k 2 k 3 , the folded shape k 1 k 2 + k 3 , etc.These two usages of "shape" are completely independent.We urge the reader to exercise caution when seeing or using this term, and we refer the reader to [55,56] for more extensive reviews.
The local shape refers to a collection of similar shapes that peaks in the squeezed limit k 1 k 2 k 3 .A typical example of local shape is

S(ρ)
Figure 1: The shape function S of the 3-point correlator as a function of momentum ratio = k 1 /k 3 .The left and right panels show schematically the shape of the tree-level process and the 1-loop process, respectively.The solid blue, dashed blue, and dotted black curves show the full 3-point function, the signal part, and the background part, respectively.This shape is called local because it can arise from a nonlinear field redefinition ζ(x) → ζ(x) + λζ 2 (x) at late universe after inflation, where ζ(x) is originally a Gaussian random field and λ is a number.Therefore we stress that the name "local" refers to a local redefinition in the late universe rather than during inflation.In other words, this is a field redefinition local to the future boundary of inflation spacetime rather than local in the bulk.A local effect on the boundary is actually very nonlocal from the bulk perspective.
The equilateral shape refers to a collection of shapes that peak at the equilateral triangle k 1 = k 2 = k 3 .A typical equilateral shape could be This shape is often generated by a local interaction in the bulk, such as (ϕ ) 3 .(Here a prime denotes conformal time derivative d/dτ ; see below.)Therefore, in our parametrization of the shape function in the squeezed limit (2), a local shape will have N = −1 while an equilateral shape has N = +1.It is clear that the local shape is more prominent in the squeezed limit.There is of course no one-to-one correspondence between the interaction types and the values of N .For example, both (ϕ ) 3 and ϕ (∂ i ϕ) 2 would give the equilateral shape N = 1.But this number still gives useful information about the interaction of the external modes.
Nonanalytic Piece.Although the purpose of this paper is to calculate the full three-point correlator for arbitrary kinematics, we will nevertheless pay special attention to the nonanalytic piece which encodes the on-shell particle production.This part is the main focus of cosmological collider physics, and for this reason, we will sometimes call it the "signal" while calling the analytic part the "background." As mentioned above, the nonanalytic piece can be characterized by a 4-parameter set (B, L, ω, ϕ).In principle, all these parameters can be calculated for specific processes and can be measured.It turns out that the scaling parameter L and the oscillation frequency ω have simple Table 1: The parameters of nonanalytic piece in (2) in several representative processes mediated by particles with spin s, mass m, and chemical potential µ.The value of B only includes leading dependence on m and µ in the large mass limit.In this table, we set H = 1.

B
L ω parametric dependence on the model parameters such as particle's mass and spin.The overall amplitude B and the phase ϕ, on the other hand, may have complicated dependences on these parameters.However, in some parameter regions, it is easy to estimate the leading dependence of B on parameters such as mass and chemical potential.We summarize some known examples of these parameters in Table 1.This table is far from exhaustive.We use it only to illustrate how the parameters of intermediate particles (spin s, mass m, chemical potential µ) are related to observables.In particular, we see that it is in general possible to tell the difference between the tree-induced and loop-induced signals, by looking at the value of L. This is essential because the value of L measures how fast the intermediate particles are diluted by inflation.In the loop process, the particles are produced in pairs and are thus diluted faster.By examining the dilution rate of intermediate modes and how they enter the signal, we can derive, for a typical tree-level or 1-loop exchange of massive particles, We plot schematically the expected shape functions from a typical tree-level process and a typical 1-loop process in Fig. 1.In this figure, we take the "background part" to be identical to the equilateral shape function (4), while the "signal part" takes the form of (5).The overall and relative amplitudes of the background and the signal are taken to be arbitrary in Fig. 1, although in given models they should be fixed and loosely related to each other.

The Formalism
We now briefly introduce the formalism that will be used for our numerical calculation in the following sections.The goal is to calculate the n-point correlators of the curvature perturbation ζ, or equivalently, the correlators of the inflaton perturbation ϕ.Upon quantization, this correlator can be interpreted as the quantum expectation value of operator products ϕ k 1 (τ ) • • • ϕ kn (τ ) at the end of inflation.Throughout the paper, we work with the dS metric ds 2 = a 2 (τ )(−dτ 2 + dx 2 ) expressed in the conformal coordinates (τ, x), and the scale factor a(τ ) = 1/(−Hτ ).In particular, the conformal time τ → (−∞, 0) and thus the end of inflation can be thought of as the τ → 0 limit.We assume the standard Bunch-Davies condition for the initial state |BD .So the npoint correlator can be written as , where τ f → 0 is the conformal time of the future infinity.Given a field theory model, this correlator can be calculated using the well-known Schwinger-Keldysh formalism.See [18] for a review and here we only summarize the main ingredients essential to our calculation.The key observation is that the expectation value can be viewed as an "in-in" amplitude, and thus can be recast into a product of two "in-out" amplitudes with the out-state scanning over a complete basis of the Hilbert space.Each of the two in-out amplitudes is amendable to a familiar path integral representation.So the correlator can be expressed as a path integral over two sets of field variables, one goes forward in time (denoted with a '+' index) and the other backwards in time (denoted with a '−' index): As above, we have two sets of fields ϕ ± , with identical action S[ϕ ± ], but with an additional minus sign in front of S[ϕ − ] to account for the "backward time".The two sets of fields are demanded equal at the future infinity by the δ-function at τ = τ f , as a consequence of summing over the out state.
The procedure then is very similar to the usual Feynman diagram expansion of the path integral, with only a few differences which we summarize now: First, each interaction vertex in a Feynman diagram is labeled by an SK index a = ±, corresponding to the two field variables ϕ ± .The minus type coupling has an additional minus sign in the vertex, coming from the minus sign in front of S[ϕ − ] in the path integral.
Second, each propagator is labeled by two SK indices at the two ends, denoted as G ab , and therefore we have 4 types of "bulk" propagators.On the other hand, if one endpoint of a propagator sits at the future boundary, then due to the identification ϕ + (τ f ) = ϕ − (τ f ), the SK index at this boundary vertex does not matter.So we will have only two types of "bulk-to-boundary" propagators, denoted by G a .
Third, it is convenient to go to the 3-momentum space thanks to the 3-dim rotation and translation symmetries.However, we do not Fourier transform the time direction, so this leads to a "mixed" version of Feynman rules.For example, each interaction vertex is associated with a 3-momentum conservation δ-function, together with an integral over time τ .Similarly, any loop in the diagram is associated with a 3-momentum loop integral, rather than a 4-momentum loop integral.
The rest of the diagrammatic rules are pretty similar to the usual Feynman rules.We refer the readers to [18] for more discussions on various technical details.

Bosonic 1-loop Process at the Cosmological Collider
The main goal of this paper is to calculate the cosmic correlator at 1-loop, with a special focus on the oscillatory signal from the loop.This section is thus devoted to a discussion of models and the related 1-loop process.This provides not only the physical motivations for our study but also the specific amplitudes that we are going to calculate in the next section.
As previous studies showed (see e.g.[9]), it is nontrivial to produce large oscillatory signals at the cosmological collider.The difficulty lies in the fact that large inflaton-matter couplings often render the matter particles too heavy to be produced, while smaller couplings reduce the signal due to vertex suppression.For many effective couplings, there is no viable parameter space in between.
The problem of getting large signals is more acute for loop processes, not only because of the additional loop factor 1/(4π) 2 .In a minimal scenario, the production is supported by the inflationary expansion and is efficient for particles with mass m H.When m H, the production is suppressed by a Boltzmann factor e −πm/H .This could introduce a suppression to the CC signal.For tree-level processes, this suppression might be tolerable and a visible signal could still be produced for m not much larger than H.But for loop processes, simple analytical estimates show that the signal is doubly suppressed, namely, by a factor of e −2πm/H .This would in general make the loop signal too small to be seen.Therefore it would be desirable to consider other production mechanisms beyond the minimal scenario.
In a well-motivated class of models, the rolling of the inflaton can produce massive particles more efficiently.See [9] for discussions.In such models, the production rate is controlled by the rolling of the inflaton through a dimension-1 parameter µ ≡ φ0 /Λ, where Λ is a cut-off scale and the perturbativity requires Λ > φ1/2 0 .In typical inflation models φ1/2 0 60H so that the production scale µ can be as high as 60H.For comparison, in the minimal scenario, the particle production is controlled by the cosmic expansion and is thus around the scale of Hubble H.
The production of massive particles via inflaton rolling can be naturally realized for particles with non-zero spin.In such cases, axion-like couplings between the inflaton and the massive fields naturally lead to particle production at the scale of µ.This includes the dim-5 couplings to a fermion (∂ µ φ)Ψ † σ µ Ψ and to a gauge boson φF F .When evaluated with a rolling inflaton background, such operators become the number density of the corresponding matter fields weighted by the helicity of the state, with the coefficient acting as a kind of "chemical potential."The size of this chemical potential is µ = φ0 /Λ and this explains why we have a new particle production mechanism at the scale of µ.The chemical potential enhanced particle production works only for one transverse polarization state so the corresponding spin-1/2 and spin-1 particles have to be produced in pairs.Thus their CC signals appear first at 1-loop order.(It is possible to have tree-level signals from longitudinal gauge bosons but it receives no enhancement from chemical potential.)Analytical estimates from previous studies have shown that these signals could be potentially large enough to be observed in the near future [10].However, a complete calculation of such a 1-loop process was not known.In this paper, we shall focus on the 1-loop process of the gauge boson and the result of the 1-loop fermion diagram will be presented in future work.
There is another class of models where the loop process could generate large signals.In such models, the curvature perturbations are generated by an additional source other than the inflaton fluctuation.Known examples of this sort include the modulated reheating scenario and the curvaton scenario [21,25].The previously mentioned problem of double suppression is partially compensated by introducing stronger coupling between the massive particle and the curvature perturbation.This cannot be realized in minimal slow-roll inflation because such a strong coupling would be inconsistent with the perturbativity.

Scalar Loop
To set the stage, we first consider the minimal case where the loop process is dS covariant.As we shall see below, the oscillatory signal in this process is in general very small.But this case is technically simpler than the more interesting cases.And also, the process is dS covariant and thus may be of theoretical interest.We are not aware of any complete analytical or numerical results for this process.The closest analytical result we have seen in the literature is a computation of 1-loop correction to the propagator in Euclidean dS completed in [57].It seems not quite straightforward to implement their result for a numerical calculation in real-time dS.So a brute-force computation directly done in real-time dS could still be useful.
We will also consider the case with chemical potential which is observationally most interesting.However, the full dS symmetry is lost in this case, so the known analytical techniques would not apply.And it seems that our numerical approach is most appropriate in this case.
We introduce a general massive scalar particle σ with bare mass m 0 , which couples to the inflaton through the dim-6 operator ( Here V (φ) is the inflaton potential that generates a slow roll motion φ1/2 0 60H const at the background level.This rolling generates a mass correction to σ so that the σ has an effective mass m2 = m 2 0 + φ2 0 /Λ 2 .To avoid the Boltzmann suppression (see below) we require m H.When both terms in the effective mass are positive, this implies that both terms should be at most of The resulting signal from this choice of parameter is thus very small since the 1-loop process is suppressed by (H/Λ) 4 .
For the current study, we will ignore the issue of this coupling suppression, insisting that the effective mass m is of O(H).It is however possible to cook up models free from both Boltzmann suppression and coupling suppression.One example is the modulated reheating as mentioned before.Another possibility is that the dim-6 operator has a flipped sign than it is in (7).Then it is possible to have both m 0 and φ2 This of course involves fine-tuning at the EFT level.But one can imagine that the theory contains a dense spectrum of scalar particles with masses around φ0 /Λ and mass differences being O(H) or smaller.Then the tuning is automatically implemented.One can further introduce a quartic self-interaction for these massive scalars to bound those states with negative effective mass.We leave a concrete model building of this sort to future works.
It is clear that the leading order contribution to the 3-point inflaton correlator from σ field is at 1-loop order, The diagram is ultraviolet (UV) divergent and needs to be regularized, which we will discuss below.
In the numerical calculation in the next section, we will also include a nonzero chemical potential to the scalar field.In realistic model building, one needs additional symmetry breaking in order to introduce a chemical potential to the scalar field at the level of mode functions.For example, [58] considered a possibility with a background static electric field acting as a chemical potential to the scalar field.This setup breaks rotation invariant and the resulting chemical potential is anisotropic.We are currently not aware of any scalar potential that does not break additional spacetime symmetries and works at the level of mode functions 3 .So an isotropic chemical potential of the scalar field can only be viewed as a toy example.

Gauge Boson Loop
The second process we are going to consider is the 3-point correlator mediated by a massive gauge boson at 1-loop.The gauge boson obtains its mass from a Higgs.We find it technically simpler to couple the gauge boson to the inflaton indirectly through the Higgs.Therefore, we consider the following Lagrangian That is, we consider a scalar QED with U (1) scalar Σ = (σ + iπ)/ √ 2 that couples to the inflaton φ through the dim-6 operator introduced in the previous subsection.The difference from the previous case is that the Higgs field Σ now has a nonzero vacuum expectation value v 2 ≡ σ 2 = φ2 0 /(λΛ 2 ).Consequently, there will be a 2-point mixing between σ and the inflaton fluctuation ϕ.This is essentially the model considered in [10] except that we do not include the dim-5 operator φF F here.We will include this dim-5 operator later.One can also include a bare mass term for Σ which in general does not alter the qualitative picture.We neglect this term here for simplicity and refer the readers to [10] for a more complete treatment.
Here a technical remark needs to be made.As long as the oscillatory signals are the only concern, we can work in the unitary gauge.The oscillatory signals are from on-shell particle production and the corresponding loop integral is never UV divergent.More explicitly, we know that the UV divergent pieces in our models can always be subtracted by local counterterms and these local counterterms always generate analytic dependence on external momenta.Therefore, the nonanalytic momentum dependence in the signal must be free of UV divergence.
In our current work, however, we are aiming at a complete calculation of the bispectrum, including both the smooth "background" and the oscillatory signal.So we need to pay attention to the UV part since the background part is UV divergent just like its flat-space counterpart.Then the calculation would be tricky in the unitary gauge.Therefore we need to consider a general R ξ gauge condition G = ∇ µ A µ − ξgvσ where g is the gauge coupling and ξ is a gauge parameter.After evaluating the Lagrangian in the dS background and standard quantization procedure, we find the Lagrangian for the fluctuating fields as where c is the ghost field.We have taken the Feynman-'t Hooft gauge ξ = 1 and neglected interaction terms irrelevant to us.All scale factors have been spelled out explicitly and so the spacetime time indices are raised or lowered by Minkowski metric η µν , except in the kinetic term of the gauge field, where = g µν ∇ µ ∇ ν .
Then at the 1-loop level, we can consider the gauge boson signal from the following two diagrams.
Here we only show the gauge boson diagrams.There are corresponding diagrams with π and ghost fields.We can drop some of these diagrams by considering special limits.For example, when λ g 2 , we only need to include π-graphs.But this would be identical to a scalar loop considered in the previous subsection and nothing new could be learned.We can also consider a different limit with λ g 2 and v H.In this limit, we can drop all graphs containing π, and retain diagrams with the gauge boson loop and the ghost loop.Of course, there can be more interesting cases if the gauge boson is coupled to external line differently (i.e., through derivative couplings).We leave those possibilities for future study.
There are additional complications in the two diagrams in (11) which are usually nonessential to the CC signals but increase the difficulty of numerical calculation.One is the presence of the dashed lines in the external legs in all diagrams.We can approximate these dashed lines by their EFT limit 1/m 2 , assuming its mass is greater than H. Another complication is the triangle loop in the second diagram in (11), which requires one more layer of time integral.In this work we shall focus only on the left diagram in (11) which we call the pinched diagram, leaving a complete evaluation of the right diagram (the triangular diagram) for future work.Incidentally, if we are only concerned with the signal part of the result, then we can approximate one of the three internal lines in the triangular diagram by its EFT limit.Thus generated "pinched coupling" is a result of integrating out one of the three internal lines, and we expect that the main contribution to this integral comes from the region when the two endpoints are separated by a distance within Hubble radius.The resulting graph is then identical to the left (pinched) diagram in (11).
Finally, we add the following operator into the Lagrangian (9) which acts as a chemical potential for the helicity of gauge bosons, This additional term modifies the dispersion relation of the gauge boson modes during inflation and thus can enhance the signal.The relevant diagrams are still the two in (11); only the propagators of the loop lines need to be replaced by the one including the chemical potential µ = φ0 /Λ F .We refer readers to [10] for more details.

Numerical Implementation and Result
In this section, we present a detailed treatment of the numerical integration.Although we focus on the pinched diagram throughout the paper, the procedures and techniques introduced here can be generalized to other 1-loop computations.We first present the expressions for the 1loop integrals and then discuss their numerical implementations.Finally, we present the numerical results.

One-Loop Integral
All the diagrams from the last subsection, including the triangle diagram after pinched approximation, can be put into the following form where the last line represents momentum permutations, which appear because of the three possible ways to connect the loop with the three external legs.In Fig. 2 we show the diagram corresponding to the first line of ( 13).Here we have suppressed all the coupling dependence that is nonessential to numerical computation.The related CC signal is given by As explained above, we need to assign SK indices a and b to the two vertices, respectively.The factor |Hτ 1 | α and |Hτ 2 | β in (13) come from the √ −g factor in the Lagrangian < l a t e x i t s h a 1 _ b a s e 6 4 = "  q k 3 < l a t e x i t s h a 1 _ b a s e 6 4 = " I h m O S 5 Z S w L Q W 8 j I 8 X v c + N E 0 l t H c = " > A A A B 7 X i c b V B N S w M x E J 3 U r 1 q / q h 6 9 B I v g q e y K q M e i F 4 8 V 7 A e 0 S 8 m m 2 T Y 2 m 6 x J V i h L / 4 M X D 4 p 4 9 f 9 4 8 9 z m E P 0 C f P 8 + s j 0 w = < / l a t e x i t > q < l a t e x i t s h a 1 _ b a s e 6 4 = " v 8 M + 7 F h s 7 r j i 4 f 8  as well as additional g µν factors in the vertices.We use G to denote the propagator of the external (inflaton) lines: In all cases, we include the derivative coupling to the inflaton by using ∂ τ G a (k; τ ) for the three external legs.For the scalar model, there is in principle a contribution from spatial derivative in the form of ∂ i G∂ i G at τ 1 -vertex.We will drop this term for simplicity although the contribution from ∂ i G∂ i G is of the same order as ∂ τ G∂ τ G.The loop integrals I ab are given by where D ab (k; τ 1 , τ 2 ) are the bulk-to-bulk propagators with where u(τ, k) is the mode function for the loop particle, either scalar or gauge boson.It turns out that the four loop integrals (I ++ , I +− , I −+ , and I −− ) are all related to the following unique integral by including various time ordering and complex conjugation: The detailed form of D > will be given below.We will compute numerically the following cases: Scalar loop.As mentioned above we will include a fictitious chemical potential µ to enhance the signal.In this case, we have α = 2, β = 3 where α, β are defined in (13), and where W is the Whittaker W function, µ ≡ µ/H and ν ≡ (m/H) 2 − 9/4.We only consider scalars with mass m > 3H/2 so ν is always real.
Gauge boson loop without chemical potential.For this example, we will stay in the Feynman-'t Hooft gauge (ξ = 1).We only consider the limit λ g 2 and v H, so that the Goldstone loops decouple, and we only need to include the gauge boson loop and the ghost loop.It turns out that the ghost loop and the A 0 loop contribute the same fictitious oscillation pattern at superhorizon scales and the two contributions cancel each other.So only the three physical degrees of freedom in the massive gauge boson contribute to the signal.In this case, we have α = 0, β = 1.The propagator will be taken to be, iν (−kτ 1 )H (2) where H (1) (H (2) ) is the Hankel function of the first (second) kind.For the gauge boson of spin-1, we have ν = (m/H) 2 − 1/4.The indices of propagators should thus be contracted in the integrand of the loop integral (17), namely, we replace D > (q; τ 1 , τ 2 )D > (p; τ 1 , τ 2 ) → D µν> (q; τ 1 , τ 2 )D µν > (p; τ 1 , τ 2 ).Again the indices here are raised by flat Minkowski metric.The propagator being proportional to g µν is a consequence of choosing ξ = 1 gauge.This naively introduces g µν g µν = 4 propagating degrees when contracting the loop propagators.One of these 4 propagating degrees is nevertheless subtracted by the ghost loop.So we will further take Gauge boson loop with chemical potential.In this case, we still have α = 0, β = 1.Due to the presence of a nonzero chemical potential, one transverse polarization will be exponentially enhanced relative to the other two polarizations.Therefore we can keep only this enhanced helicity state to a good approximation.The propagator in this case is where e (h) (k) is the polarization vector for helicity state h.
In all cases above, the loop integral is divergent in the UV just like in flat space.However, the asymmetric treatment of space and time makes it difficult to isolate and extract the divergence directly from the integrand.Therefore we have to choose a regularization method that is straightforward for numerical implementation.We will adopt a Pauli-Villars regulator, by rewriting the loop integral (17) as where the regularized propagator D reg is the original propagator with mass m subtracted by the same propagator but with a heavier mass M , Then the heavy mass M serves as an effective cutoff of the loop integral.
After the regularization, one still needs a renormalization scheme to fix the answer of the loop correction.This can in principle be done in given models.Since this part is less relevant to our numerical computation and is more model-dependent, we will only comment briefly on it here.Take the scalar loop model (7) as an example, one needs to include a counterterm in order to cancel the regulator dependence in the loop integral.Then a renormalization scheme is needed to fix the regulator-independent part of the counterterm coefficient δ C .As we will see below, the nonanalytic signal part of the correlator is free from the UV divergence and is independent of the renormalization scheme, since it is essentially from the on-shell particle production.
Therefore it is at least in principle possible to extract the mass and the coupling of the intermediate particles by measuring the oscillatory signal.On the other hand, the background part of the signal receives UV contributions and its loop correction will be renormalization-dependent.The renormalization condition can thus be determined by measuring the overall amplitude of the background part of the correlator, namely the coefficient A in (2).Then, using the measured mass, coupling, and amplitude A, one can fix the coefficient δ C .To tell apart the background from the signal, one needs to measure a range of external momentum configurations, namely, to measure the shape dependence.After the above subtraction is done, one can then use the loop diagram result plus the counterterm to predict the shape function for wider range of shapes.While this procedure is in line with the usual renormalization story in flat space, an important difference needs to be noted: In flat-space QFT, changing the external momenta amounts to changing the energy scale.For inflation correlators, however, the external momenta only label the comoving scales but not the physical energy scales.Therefore, changing external comoving momenta only probes the shape dependence in a scale invariant model, whereas the process always happens at a fixed physical scale (namely the scale of particle production at either the chemical potential µ or Hubble scale H.) A complete set of renormalization conditions of course requires more careful treatment of various loop diagrams contributing to inflaton's correlators.The above simple argument nevertheless tells us that the background part of the 3-point function, in general, cannot be calculated unambiguously without a careful renormalization procedure.On the other hand, the signal part is free from UV issues.Below we will see that this is indeed the case by showing that the signal part is regulator-independent while the background part is regulator-dependent.Therefore, the background part of our calculation below should be treated only as an indication of the overall size of the regularized loop diagram.One should not view them as the direct prediction of the model and should not compare them with data directly.

Setup and Procedure
We separate the numerical evaluation of the integral (13) into two stages: (i) For each combination of SK indices (ab), we generate a discrete 2-dim grid of the conformal time (τ 1 , τ 2 ) and compute the loop integral ab on each grid point (with Wick rotation as introduced below).We parallelize the evaluation of the grids on a cluster.This stage is computationally expensive, given the grids are tightly spaced.We utilize several properties to reduce the computation load.First, given the reality condition of SK diagrams, we only need to compute the grids with SK indices of (ab) = (++) and (+−).No evaluation for (ab) = (−−) or (−+) is needed.Second, as explained in detail later, we only need to evaluate the τ 1 ≥ τ 2 part of the (ab) = (++) and (+−) grids, which reduces about half of the computation task.
(ii) After obtaining the tabulated grids of (τ 1 , τ 2 , I ab ) with (ab) = (++) and (+−) from Stage (i), we dress ab with the bulk-to-boundary propagators and other pre-factors as shown in (13).We then interpolate the dressed grids and integral over τ 1 and τ 2 to get the 1loop result, B ab (k 1 , k 2 , k 3 ), for a particular pair of SK indices (ab) and external momentum configuration (k 1 , k 2 , k 3 ).Since the time grids are often tightly spaced, we can approximate the time integral by a Riemann sum of the dressed grids to speed up the evaluation.The final result is given by summing over all combinations of SK indices as well as the momentum permutations, i.e., B(k , where we applied the reality condition in the second equality.This stage does not require many computational resources and can be finished in a much shorter period compared to Stage (i).
For the scalar or gauge boson loop we considered, the loop integrals in Stage (i) contain products of Whittaker functions (e.g.Eqs. ( 18) and ( 19)) with complex orders, i.e., W κ,λ (z) with κ, λ ∈ C. It is convenient to implement the numerical code in Mathematica 12, which has native support for those functions.An alternative choice is to use mpmath [59], a Python library.In App.A, we crosscheck results from Mathematica 12 code with those from mpmath 1.1 code for the CC signal of the tree-level process for complex scalar in an electric field.We find a good agreement between the two.Further relevant details for the numerical evaluation are listed below: • The choice of chemical potential and mass parameter.The Whittaker functions, W ∓iμ,±iν (±2ikτ ), are highly oscillatory when µ or ν is large.As a result, numerical evaluations for the loop integral ab , in general, takes a longer time as μ or ν increases.Here we choose μ and ν to be O(1).For the Pauli-Villars regulator (22), we set the mass parameter M such that νreg of the Whittaker functions W ∓iμ,±iνreg (±2ikτ ) inside D M (k, τ 1 , τ 2 ), satisfies νreg ≥ 2ν.
• Range of the magnitude of the external momentum.We will present the momentum dependence of our results in two ways.One is the squeezed limit (k , where the CC signal is expected to appear.The other is a more general range of momenta, the near-equilateral limit (k , which is often adopted in literature for reporting the bispectrum.Signals with k 1 = k 2 ≤ k 3 in the near-equilateral limit can be viewed as a natural extension of those with k 1 = k 2 k 3 in the squeezed limit.By the scale invariance, we can fix k 3 = 1 without loss of generality.To set the range of the time grid for Stage (i), we first need to specify the range of k 1,2,3 of interest.For the squeezed limit, we consider k ∈ [1, k max ] with k max ∼ 10 3 , where we expect to see O(ν) oscillations per decade The loop integral values for the grid with τ 1 < τ 2 (white part) can be copied from the evaluated grid up to the exchange of τ 1 ↔ τ 2 .A similar reduction of computational task can be achieved for the (ab) = (+−) diagram, where we only need to evaluate the grid with τ 1 ≥ τ 2 .See text for more details.(Right) An overlay of the integration contour of τ on top of the complex plot for the Whittaker function W −4i,4i (2iaτ ) with SK index a = +.The color codes the features of the complex function [60] and the white wedge along the positive imaginary axis represents the branch cut of the Whittaker function.To achieve a better convergence for the time integral while avoiding the branch cut, we rotate τ where is a small positive real number.in k.For the near-equilateral limit, we consider the triangle region with k 1,2 ∈ [k min , 1] with k min ∼ 0.1.
• Range and spacing of the time grid.In (13), we take τ 1,2 ∈ (−∞, 0) for the time integral.But for the numerical evaluation, the integral range should be finite.Therefore we need to introduce a finite interval To capture all relevant physics, the initial time τ i needs to be set early enough so that all relevant modes are deep inside the horizon (|kτ i | 1) and well before the particle production (|kτ i | μ).The final time τ f needs to be late enough such that all modes are well outside the horizon (|kτ f | 1).Given k ∈ [k min , k max ] with k min ∼ 0.1 and k max ∼ 10 3 , we choose τ i = −200 and τ f = −10 −5 and log-evenly spacing the interval τ 1,2 ∈ [τ i , τ f ] into N grid = 440 pieces.In the end, we get a 2-dim time grid of (τ 1 , τ 2 ) with (N grid + 1) 2 grid points in total for each SK indices combination.Given the setup, each grid point occupies the same area in the log space, [(ln |τ i | − ln |τ f |)/N grid ] 2 , which we will use as the measure for the Riemann sum in Stage (ii).
• Reducing the grid evaluation for the diagram with (ab) = (++).For both scalar and gauge boson loop diagrams, the key part of the CC signal is given in the form of where we only keep track terms of interest.For the (ab) = (++) diagram, the loop propagator is a piece-wised function, D ++ (q; τ 1 , τ 2 ) = D > (q; τ 1 , τ 2 )θ(τ 1 − τ 2 ) + D < (q; τ 2 , τ 1 )θ(τ 2 − τ 1 ), as shown in Eq. ( 16).Using the relation D < (q; τ 1 , τ 2 ) = D > (q; τ 2 , τ 1 ), the integral ( 24) can be expressed as where the two terms inside the square bracket of the last line of ( 26) are identical with the exchange of τ 1 ↔ τ 2 .The same symmetry can be found when replacing the propagator D ++ with the regulated propagator D reg,++ .Given the symmetry between the two terms, we only need to evaluate I After evaluating the loop integrals and obtain the grid (τ 1 , τ 2 , I ++ ) with τ 1 > τ 2 , one can simply switch τ 1 and τ 2 and copy the loop integral value I ++ to tabulate the rest of the grid. 5We illustrate the procedure in the left panel of Fig. 3.
One subtlety arises for the parameter space just around the diagonal axes, τ 1 = τ 2 .It is the place where the oscillations in the propagators are significant.If we strictly select grid points that satisfy τ 1 > τ 2 , the parameter space will be left out given the finite spacing of the time grid.In the code, we include grid points with τ 1 = τ 2 , illustrated as darker colored squares in the left panel of Fig. 3, to represent the values of I ++ around them. 6   • Change variables for the time integral.Because the time grids are log-evenly spaced, it is natural to consider s 1,2 ≡ ln |τ 1,2 | as the integration variables for the time integrals in Stage (ii).Such a change of variable is convenient for the Riemann sum approach given the measure for each grid point is a fixed value, ∆s The change of time variables is also helpful for the interpolate-then-integral approach because it improves the sampling of the integrand for the time integral.
5 Note that although I +− (k 3 , τ 1 , τ 2 ) shares the same integrand as I ++ (k 3 , τ 1 , τ 2 ) with τ 1 < τ 2 , we could not simply copy the values between the two given their differences in the Wick rotation treatment: we rotate τ 1,2 → −iτ 1,2 (τ 1 → +iτ 1 and τ 2 → −iτ 2 ) when evaluating I +− ).See below for Wick rotation. 6Incidentally, we note that the τ 1 τ 2 region contributes most to the "background" of the correlator and contributes little to the "signal."Therefore one can imagine to approximate the contribution from this region by taking another pinched limit (namely the EFT limit).While this is an interesting point to check numerically in a future work, we must note immediately that this statement is not true in general.For the 4-point correlator with s-channel exchange, for example, the τ 1 = τ 2 region can contribute significantly to the signal as well, depending on the choice of external momenta.
• Wick rotation of the time variables.The time integral in Stage (ii) contains a lot of oscillations from the bulk-to-boundary propagators ∂ τ G ± (k, τ ) ∼ e ±ikτ , which is numerically difficult to handle.Hence we perform the Wick rotation on τ 1,2 , for computations in Stage (i) and (ii), to improve numerical convergence in the early-time limit (τ → τ i ). 7Thanks to the regularization, our loop integrand vanishes in the large loop-momentum limit, as shown in App.B. We need to be careful to avoid branch cuts in the Whittaker functions as described below.The sign of the rotation, whether τ → +iτ or τ → −iτ , is determined by the i -prescription, which is practically equivalent to such that the exponential in the bulk-to-boundary propagators are suppressed when |τ | → ∞, i.e., e ±ikτ → e −k|τ | = e +kτ .Therefore we take τ → −iτ (τ → +iτ ) if τ is companied with SK index a = + (a = −).
• Avoiding the branch cut.The Whittaker W function has a branch cut emanating from 0 to ∞.In the implementation of Mathematica 12 and mpmath 1.1, the branch cut of the Whittaker function W κ,λ (z) is along the negative real axis of z, which should be avoided when performing Wick rotation.This yields a problem for the Whittaker functions with the argument z = 2iakτ when the SK index a = +.After performing the Wick rotation described earlier, z = 2iakτ = 2ikτ → z = 2kτ ∈ R − , which hits the branch cut of the Whittaker function.To avoid the branch cut, we choose to downshift the wick-rotated z by a small imaginary number 2ik with ∈ R + such that z = 2ikτ → z = 2kτ − 2ik (or equivalently τ → −iτ − ).The right panel of Fig. 3 is an illustration of the rotation procedure for W −4i,4i (2iaτ ) with a = +.Given the range of momentum and time we considered, we fix = 10 −9 and perform the extra shift after rotating τ 1,2 , i.e., τ 1,2 → −iτ 1,2 − , when evaluating ++ . 8  • Reducing the grid evaluation for the diagram with (ab) = (+−).Under complex conjugation, the Whittaker function has the following property We utilize it to reduce the computational load for the (ab) = (+−) diagram.The CC signal for the diagram is given by where we rotated the time variables according to the Wick rotation rules described above. 7There is a potential IR problem of full Wick rotation for processes involving mutual cancellation of IR divergences among different diagrams.Our diagrams are free from such problems. 8For I +− , the Wick-rotated τ 1,2 do not hit the branch cut under the rotation rule.Hence the extra small shift is not necessary.
For each propagator, we found where • • • represents pre-factors that are symmetric under τ 1 ↔ τ 2 .We apply Eq. ( 27) in the second line of Eq. ( 29).The same conjugation property is valid for D reg,> .Using Eq. ( 29), the signal can be split into two "symmetric" parts where the two terms inside the curly bracket of the integral (30) are identical up to the exchange of τ 1 ↔ τ 2 and the complex conjugation.The same symmetry can be found when replacing D > with D reg,> .As a result, similar to the case for the (ab) = (++) diagram, we only need to evaluate τ 1 ≥ τ 2 part of the time grid for the loop integral I +− .After evaluating I +− for τ 1 ≥ τ 2 , one can simply switch the position of τ 1 and τ 2 and conjugate the corresponding I +− values to tabulate the rest of the time grid.The first term in the loop integral (26) and that of (30) can be combined together as a unified function I (1) for (ab) = (++) and t 1 = −i|τ 2 | and t 2 = i|τ 1 | for (ab) = (+−).For each grid, we need to evaluate (N grid + 1)(N grid + 2)/2 grid points.Together (N grid + 1)(N grid + 2) grid points are needed to be evaluated in order to get the CC signal for a pinched loop diagram.This is about half of the original total number of grid points, 2(N grid + 1) 2 .
• Integral of the loop momentum.For the pinched loop diagram, the momentum running in the loop propagators are q and p ≡ q − k 3 respectively.We have the freedom to build the coordinate for q and we choose a spherical coordinate with the polar direction (θ q = 0) aligning with the direction of k 3 .Under such a coordinate, the magnitude of p is given by p = q 2 − 2qk 3 cos θ q + k 2 3 .And the integration over q can be expressed in spherical coordinates d 3 q/(2π) 3 → +1 −1 d cos θ q +∞ 0 dq q 2 /(2π) 2 , where the azimuth angle of q has been integrated out.
• Patched Whittaker function.It is quite time-consuming to numerically integrate the loop integrand over q ∈ [0, ∞).To speed up the evaluation, we adopt an approximation for the Whittaker functions, dubbed as the patched Whittaker functions W κ,λ (z): we replace the Whittaker functions W κ,λ (z) with its asymptotic expansion at z = 0 and z = ∞ for |z| |λ| and |z| |λ| respectively, while keeping the full function for the intermediate range.To be concrete, we use the third (second) order expansion in z for the region |z| |λ| where and In App.A, we crosscheck the CC signals from the tree-level diagram for complex scalar in an electric field using the patched Whittaker functions (31) with those using the full Whittaker functions.We find a good agreement between the full numerical Whittaker and the triply patched Whittaker.Besides we also find a good agreement between the CC signals computed in both ways.The numerical results also agree well with the results from the analytical expression for the (ab) = (+−) diagram.
• Contributions from the momentum permutations.In the above computation, we focus on the case where the loop integral I ab is controlled by a fixed external momentum k 3 = 1.To get the contribution from the momentum configuration (k 1 ↔ k 3 ) or (k 2 ↔ k 3 ), we do not need to re-evaluate the loop integral controlled by k 1 or k 2 .Instead we can utilize the scale invariance of ).The contributions from the two momentum permutations, S ab (k 1 ↔ k 3 ) and S ab (k 2 ↔ k 3 ), are respectively given by and where we rescale k 1 , k 2 , and k 3 of S ab (k 1 ↔ k 3 ) (S ab (k 2 ↔ k 3 )) by a common factor of k 3 /k 1 (k 3 /k 2 ) in the derivation of Eq. (34) (Eq.( 35)).Below we will use S ab (k 1 , k 2 , k 3 ) to represent the sum of the signal from the momentum permutations.Utilizing the above procedures, parameter choices, and special treatments, we find the loop integral at each time grid point in Stage (i) takes a few CPU hours to finish on average.Given a total number of 194,922 grid points (after applying the grid reduction) to evaluate, it costs O(10 5 ) CPU hours to get the final CC signal for one pinched loop diagram.

Results
The 1-loop 3-point correlator is the sum of contributions from all SK diagrams, S(k 1 , k 2 , k 3 ) = a,b=± S ab = 2Re S ++ +2Re S +− , where we applied the reality condition in the last equality.Given the reality of the dressing functions and pre-factors, the value of Re S ++ (Re S +− ) is determined by the real part of the loop integral I By taking the time integral (Riemann sum) over the grids of the dressed loop integral and summing over the momentum permutations, we obtain the CC signal for a given momentum configuration (k 1 , k 2 , k 3 ).In the left panel of Fig. 5, we show the resulting signal, rescaled by k 1 /k 3 , in the squeezed limit, k 1 = k 2 k 3 .The resulting signal demonstrates periodic oscillations as a function of the logarithm of the momentum ratio, log 10 (k 1 /k 3 ), with a period of T = π log 10 e/ν.We insert a set of gray bands, which are evenly separated by ∆ log 10 k 1 /k 3 = T , in the panel as a  guide to the pattern.Such oscillations can be also observed from the results of the diagrams with the same(opposite)-sign SK indices, 2Re S ++ (2Re S +− ), that contribute to the final signal.To see the oscillatory pattern more clearly, we attenuate the smooth component of the signal.In the lower sub-panel, we show the filtered results by passing the signal through a high-pass filter with a Gaussian window function and a cut-off frequency of ω c = 0.05ν/(π log 10 e).The amplitude of the wiggle is about one to two orders of magnitude smaller than that of the full signal.We apply the same high-pass filter for all the filtering procedures in the main text and crosscheck it with other filtering methods in App. C.
Next, we look into the dependence of the signal on the loop regulator, ν reg .The right panel of Fig. 5 shows the resulting signal from the diagrams with the same-sign SK indices, S ++ +S −− = 2Re S ++ , in the squeezed limit.The solid blue and dashed orange curves are results with parameter (μ, ν, νreg ) = (5,4,8) and (μ, ν, νreg ) = (5,4,12) respectively.The resulting signals from different values of ν reg share the same periodicity and phase, and a mild difference in their amplitudes.This indicates the divergence of the loop diagram is logarithmic.The two signals are almost identical after passing through the same high-pass filter.Now we switch to the gauge boson loop.The upper left panel of Fig. 6 shows the signal for the 1-loop diagram of gauge boson without chemical potential, where the parameter is set to (μ, ν, νreg ) = (0, 2, 4).The oscillatory signal is relatively tiny since the chemical potential is absent.(3.9, 4, 8).The oscillatory pattern is obvious in the filtered signal.The shape and amplitude of the (filtered) signal are similar to those of the scalar loop shown in the left panel of Fig. 5. Unlike the case without chemical potential, contributions from the diagrams with the same/oppositesign SK indices are comparable.The lower left panel of Fig. 6 shows the results for a slightly different parameter choice (μ, ν, νreg ) = (4, 4, 8), where we found the signal is generally a factor of exp(2π∆μ) ∼ 1.8 larger than that with (μ, ν, νreg ) = (3.9, 4, 8).We demonstrate the ν regdependence in the lower right panel, where the solid blue and dashed orange lines are the resulting total signals with parameter (μ, ν, νreg ) = (4, 4, 8) and (μ, ν, νreg ) = (4,4,12) respectively.Similar to the scalar scenario, the two signals share a mild difference in their amplitudes and yield almost identical oscillation patterns after passing through the same high-pass filter.In Fig. 7, we plot the 3-point correlator mediated by the gauge boson loop for more general  momentum configurations.The gauge boson loop is dominated by the negative helicity mode for large chemical potentials.We show the four sets of the parameter (μ, ν, νreg ) described earlier.
We use the dash lines to mark the boundaries of the regions with k 1,2 |τ i | max(μ, 1), 9 where the numerical results are not handicapped by the finite value of τ i during the time integral.

Comparison with Analytic Estimates
The loop amplitudes calculated above have been estimated in previous works.The estimates usually assumed a late-time expansion of the loop propagator which is not entirely valid for 3point functions.Explicitly, it was estimated that the signal part of a 1-loop process mediated by a boson with mass m ≡ m/H 1 and chemical potential µ ≡ µ/H takes the following form as a function of ≡ k where ν = m 2 − 9/4 for scalar and ν = m 2 − 1/4 for spin-1 gauge bosons.The phase ϕ can be calculated and depends on µ and ν in complicated ways.The powers y and z depend on interactions and are difficult to estimate precisely.It is known that late-time expansion of the intermediate propagators can capture the exponential dependence e 2π( µ− ν) correctly but fails to get the correct power dependence µ y m z .But given that the most sensitive dependence comes from the exponential factor, it is still useful to compare our numerical results with the analytical estimates.In Fig. 8, we superimpose the result of analytical estimates (36) with the filtered numerical results for four 1-loop diagrams we considered, setting y = z = 0, and tuning the phases ϕ to match the peaks and valleys of the analytical and numerical oscillatory signals for the region k 1 /k 3 10.We found that the overall amplitudes of the numerical results and analytical estimations are compatible.We also observe a small difference in the oscillation frequencies between the two.Our numerical results show a slightly frequency than analytical estimates ω = 2ν/ log 10 e in all examples.We are currently unaware of the origin of this small discrepancy.Partly due to this difference and partly due to the failure of analytical estimates at small momentum ratios, the two results disagree in phases at lower values of k 1 /k 3 , while they are in reasonable agreement for k 1 /k 3 20, as shown in Fig. 8.

Conclusions
We performed a systematic study of the numerical evaluation of the inflation correlators at the 1-loop level.In particular, we presented the first numerical results for a set of 3-point 1-loop processes with spin-0 or -1 particles running in the loop.Without applying crude approximations, our numerical results offer a set of precise templates for the 1-loop process after imposing appropriate renormalization conditions.Our results include both "signal" and "background" parts, that can be used in searches with observational data.While the "background" can be several orders of magnitude larger, we show that the signal can be separated out by applying high-pass filters.We have also compared our numerical results with analytic estimates adopted in the earlier literature.The amplitudes obtained by both approaches are compatible.At the same time, there is a shift in the frequency of the oscillatory part for k 1 /k 3 at a several-percent level.Further clarifying the origin of this shift could be a fruitful future direction to pursue.
As shown in this paper, the numerical implementation of the relevant 1-loop integrals turns out to be highly nontrivial.Many subtleties could appear along with the calculation, including the appropriate choices of integral range and gridding, the highly oscillatory integrand in certain parameter regions, the divergence of the integral and its regularization, the analytical structure of special functions that is relevant when performing Wick rotation, etc.We have spelled out all these subtleties in the paper, which could be useful for more extensive numerical studies in the future.We stress that our method does not rely on approximations that are often used in previous analytical studies and work usually only in several limiting momentum configurations.Our method thus applies to arbitrary momentum configurations.
It would be interesting to extend our current project to other interesting situations, such as the fermionic 1-loop process, more general triangular 1-loop processes (without taking pinched limit), and more general couplings (e.g.φF F mentioned in the paper).It is also straightforward to apply our method to 1-loop mediated trispectrum (4-point correlation functions).The results of these studies can be a useful check of previous analytical estimates and are also useful for building templates for future observations of the 3-point functions.
As the first attempt of numerical evaluation of 1-loop bispectrum, our method is essentially a direct implementation of Feynman integrals from first principles.As we have shown, the evaluation of relevant Feynman integrals is computationally heavy.It would be inefficient to apply our method directly for parameter scanning and templates building today.For these purposes, it would be desirable to look for both analytical simplifications and better numerical strategies.We leave these directions for future studies.
) where 2 F1 is the regularized hypergeometric function.Therefore, the tree-level diagram for complex scalars in an electric field provides an ideal scenario for us to crosscheck the convergence of various numerical implementations described in Sec.4.2.We first compare numerical results from those using the full Whittaker functions to those using the patched Whittaker functions.We then compare the numerical results from coding in Mathematica to those from coding in mpmath.

A.1 Whittaker function vs. patched Whittaker function
We consider the CC signal, S(k 1 , k 2 , k 3 ), in the squeezed limit k 1 = k 2 k 3 , with the electric field strength and the mass parameter (μ, ν) = (0.1, 4), (1,4), and (4, 4).For a given set of parameter, the total signal is the sum of positive and negative electric field strength parameter, i.e., S = ab S ab (μ, ν)+ ab S ab (−μ, ν) = S ±± (μ, ν)+ S ±∓ (μ, ν)+(μ → −μ) = 2Re[S ++ (μ, ν)+ S +− (μ, ν)] + (μ → −μ), where we apply the reality condition of SK diagrams in the last equality.The final signals, rescaled by k 1 /k 3 , are shown in the first row of Fig. 9, where the blue lines are the results from the numerical evaluation using the full Whittaker functions and the orange lines are the results using the patched Whittaker functions (31).We see a good agreement between the two implementations.Such agreement can be also seen at the level of component that contributes to the signal, as shown in rows 2-5 of Fig. 9.Note that for S ±∓ = 2Re S +− components (row 3 and 5), we also include results from the analytical expression of (38), which are in good agreement with the numerical results.
All the CC signals show oscillatory patterns with a periodicity of T = 2π log 10 e/ν ≈ 0.68 with respect to log 10 (k 1 /k 3 ).For each set of parameters, the dominant contributions are from the diagrams with SK indices of the same sign.As μ increases, the overall size of the signal increases, and the oscillatory patterns of the signals begin to dominate over the smooth backgrounds.

A.2 Mathematica vs. mpmath
We employ Mathematica for most of the computational tasks through the paper.Here we crosscheck the numerical results from Mathematica code to those from mpmath code.For each code implementation, we check results from computation with the full Whittaker functions as well as the patched Whittaker functions.The upper panel of Fig. 10 shows the resulting signals for the parameter (μ, ν) = (5,4).The results from Mathematica are in good agreement with those from mpmath for both the full Whittaker and the patched Whittaker functions.The middle and lower panels of Fig. 10 show such agreement exists at the individual SK diagram level.We again add the result from the analytical expression of 2Re S +− to the lower panel of Fig. 10, which agrees well with the numerical results.We also include the analytical result from Eq. (38) (dashed black) to compare with the numerical results.A set of gray bands, which are evenly separated by ∆ log 10 (k 1 /k 3 ) = 2π(log 10 e)/ν, are included for all panels as a guide to the oscillatory pattern.

B The behavior of the loop integrand at large loop momentum
For the scalar or vector loop integral I ab , its integrand can be expressed in terms of a single function L as L(k 3 ; τ A , τ B ) ≡ q 2 τ α A τ α B e 2μπ 4pq W −iμ,iν (2piτ A )W iμ,iν (−2piτ B )W −iμ,iν (2qiτ A )W iμ,iν (−2qiτ B ), (40) where c q ≡ cos θ q , p = q 2 − 2qk 3 c q + k 2 3 , and the power α = 2 (α = 0) for the scalar (gauge boson) loop.The loop integral from the (ab) = (++) diagram is given by I ( Again we can split L into the τ 1 > τ 2 part and the τ 1 < τ 2 part and only need to perform integral for one of them as discussed in Sec.4.2.
For our numerical procedure, it is important that the loop integral always yield finite results on a given time grid of (τ 1 , τ 2 ).Therefore, we need to check whether the loop integrand vanishes at the large loop momentum limit, q → ∞.Under such limit, the other momentum inside the loop integral p ≈ q since q −k 3 ≤ p ≤ q +k 3 .(We fixed k 3 = 1.)Meanwhile, the Whittaker function can be expanded as W κ,λ (z) ≈ z κ e −z/2 [1 + O(z −1 )] at z → ∞.Under such limit, Eq. ( 40 Let us first consider the large q behavior of the integrand of I given the suppression from e −2q(|τ 2 |−|τ 1 |) .For the same reasoning the second L term of the integral (41), L(k 3 ; τ 2 , τ 1 )| τ 2 >τ 1 , also vanishes at large q.Next, we consider the large q behavior of the integrand of I 4 h X P w 4 A r q c A c N a A I D B c / w C m + O c V 6 c d + d j M V p w 8 p 1 j + A P n 8 w f 4 w Z E g < / l a t e x i t > b < l a t e x i t s h a 1 _ b a s e 6 4 = "f L j Q A b o E 2 T x / t 8 9 m + I 1 v v s s h K d M = " > A A A B 8 X i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s x I U Z d F N y 4 r 2 A e 2 p d x J M 2 1 o J j M k G a E M / Q s 3 L h R x 6 9 + 4 8 2 9 M p 7 P Q 1 g O B w z n 3 k n O P H w u u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / + g f H j U 0 l G i K G v S S E S q 4 6 N m g k v W N N w I 1 o k V w 9 A Xr O 1 P b u d + + 4 k p z S P 5 Y K Y x 6 4 c 4 k j z g F I 2 V H n s h m r E O U p w N y h W 3 6 m Y g q 8 T L S Q V y N A b l r 9 4 w o k n I p K E C t e 5 6 b m z 6 K S r D q W d p 9 r V K / y e s o w g m c w j l 4 c A V 1 u I M G N I G C h G d 4 h T d H O y / O u / O x G C 0 4 + c 4 x / I H z + Q P 3 P J E f < / l a t e x i t > a < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 Z 2 K f h I S 9 L z z h O H r e 5 B R 7 Y j G v C 0 = " > A A A B + H i c b V D L S s N A F L 2 p r 1 o f j b p 0 M 1 g E N 5 Z E i 7 o s u n F Z w T 6 g D W E y n b R D J w 9 n J k I N + R I 3 L h R x 6 6 e 4 8 2 + c t F l o 6 4 H L P Z x z L 3 P n e D F n U l n W t 1 F a W V 1 b 3 y h v V r a 2 d 3 a r 5 t 5 + R 0 a J I L R N I h 6 J n o c l 5 S y k b c U U p 7 1 Y U B x 4 n H a 9 y U 3 u d x + p k C w K 7 9 U 0 p k 6 A R y H z G c F K S 6 5 Z H X h + + p C d 5 m 3 i n m e u W b P q 1g x o m d g F q U G B l m t + D Y Y R S Q I a K s K x l H 3 b i p W T Y q E Y 4 T S r D B J J Y 0 w m e E T 7 m o Y 4 o N J J Z 4 d n 6 F g r Q + R H Q l e o 0 E z 9 v Z H i Q M p p 4 O n J A K u x X P Ry 8 T + v n y j / y k l Z G C e K h m T + k J 9 w p C K U p 4 C G T F C i + F Q T T A T T t y I y x g I T p b O q 6 B D s x S 8 v k 8 5 Z 3 b 6 o N + 4 a t e Z 1 E U c Z D u E I T s C G S 2 j C L b S g D Q Q S e I Z X e D O e j B f j 3 f i Y j 5 a M Y u c A / s D 4 / A G t 0 Z M d < / l a t e x i t >

5 9 3 3 <
5 m L e u O P n M E f y B 8 / k D 8 s i P 7 A = = < / l a t e x i t > k l a t e x i t s h a 1 _ b a s e 6 4 = " E V Q S x / P p Y S m V 0 t D l 4 J S x G E + 9 b Z w = " > A A A B 7 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m k q M e i F 4 8 V 7 A e 0 o W y 2 m 3 b p Z h N 3 J 0 I J / R N e P C j i 1 b / j z X / j t s 1 B W x 8 M P N 6 b Y

2 Figure 2 :
Figure 2: The pinched diagram.a and b are SK indices.k 1,2,3 are external 3-momenta.τ 1 , τ 2 are conformal time variables.The black and blue lines represent the boundary-to-bulk propagators and the internal loop propagators respectively.

Figure 3 :
Figure 3: (Left) 2-dim time grid in (ln |τ 1 |, ln |τ 2 |) for the pinched 1-loop diagram with SK indices (ab) = (++).We only need to evaluate the loop integrals for the grid with τ 1 ≥ τ 2 (colored part).The loop integral values for the grid with τ 1 < τ 2 (white part) can be copied from the evaluated grid up to the exchange of τ 1 ↔ τ 2 .A similar reduction of computational task can be achieved for the (ab) = (+−) diagram, where we only need to evaluate the grid with τ 1 ≥ τ 2 .See text for more details.(Right) An overlay of the integration contour of τ on top of the complex plot for the Whittaker function W −4i,4i (2iaτ ) with SK index a = +.The color codes the features of the complex function [60] and the white wedge along the positive imaginary axis represents the branch cut of the Whittaker function.To achieve a better convergence for the time integral while avoiding the branch cut, we rotate τ ∈ [τ i , τ f ] to −iτ ∈ [−iτ i − , −iτ f − ], where is a small positive real number.

6 Figure 4 :
Figure 4: (Left) The real part of the loop integral I

4 (Figure 10 :
Figure10: (Upper ) The CC signal S, rescaled by k 1 /k 3 , in the squeezed limit, k 1 = k 2 k 3 , for a tree-level diagram for complex scalar in an electric field.The electric field strength and the mass parameter are set such that (μ, ν) =(5,4).We compare numerical results from Mathematica 12 and those from mpmath 1.1, together with those from the full Whittaker functions and those from the patched Whittaker functions.(Middle)S ±± component of S (S = S ±± + S ±∓ = 2Re(S ++ + S +− )).(Lower ) S ±∓ component of S. We also include the analytical result from Eq. (38) (dashed black) to compare with the numerical results.A set of gray bands, which are evenly separated by ∆ log 10 (k 1 /k 3 ) = 2π(log 10 e)/ν, are included for all panels as a guide to the oscillatory pattern.