# Precision Calculation of Inflation Correlators at One Loop

Lian-Tao Wanga111Email<EMAIL_ADDRESS>Zhong-Zhi Xianyub222Email:
<EMAIL_ADDRESS>Yi-Ming Zhongc333Email<EMAIL_ADDRESS>
_a Department of Physics, University of Chicago, Chicago, IL 60637, USA_
_b Department of Physics, Tsinghua University, Beijing 100084, China_
_c Kavli Institute for Cosmological Physics, University of Chicago, Chicago,
IL 60637, USA_

###### Abstract

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.

##  1 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 $\mathcal{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\sim 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 $\zeta$,
producing a characteristic oscillatory signal in various soft limits of
$\zeta$-correlators, including the “squeezed limit” ($k_{3}\ll k_{1}\simeq
k_{2}$ where $k_{i}=|\mathbf{k}_{i}|$) of the 3-point function (bispectrum),
and the “collapsed limit” ($|\mathbf{k}_{1}+\mathbf{k}_{2}|\ll k_{i}$,
$i=1,\cdots,4$) 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 [12, 13, 14, 15, 16, 17, 18,
8, 19, 20, 21, 22, 23, 24, 25, 26, 9, 27, 28, 29, 10, 30, 31, 11, 32]. 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 chemical-potential-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
axion-type dim-5 coupling to the fermions or gauge bosons [36, 37, 38, 39, 40,
8]. 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 [41, 18]. 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 symmetry-based 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 $\phi F\mkern
2.0mu\widetilde{\mkern-2.0muF\mkern-2.0mu}\mkern 2.0mu$ 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 1-loop 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.

##  2 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.

###  2.1 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 $\mathbf{k}_{i}~{}(i=1,2,3)$ and it depends on
$\mathbf{k}_{i}$ only through their magnitudes, $k_{i}\equiv|\mathbf{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}\ll k_{1}\simeq k_{2}$.

For the 3-point correlator, it is conventional to define a dimensionless shape
function $\mathcal{S}(k_{1},k_{2},k_{3})$ in the following way,

$\langle\zeta_{\mathbf{k}_{1}}\zeta_{\mathbf{k}_{2}}\zeta_{\mathbf{k}_{3}}\rangle^{\prime}=\displaystyle\frac{\,{(2\pi)^{4}\mathcal{P}_{\zeta}^{2}}\,}{k_{1}^{2}k_{2}^{2}k_{3}^{2}}\mathcal{S}(k_{1},k_{2},k_{3}),$
(1)

where a prime in $\langle\cdots\rangle^{\prime}$ means to strip away the
momentum-conserving factor
$(2\pi)^{3}\delta^{(3)}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})$, and
$\mathcal{P}_{\zeta}\simeq 2\times 10^{-9}$ is the nearly scale-invariant
power spectrum, which measures (the square of) the size of the curvature
perturbation $\zeta$. In the single-field slow-roll models, $\zeta$ is related
to the inflaton fluctuation $\varphi$ via
$\zeta_{\mathbf{k}}\simeq-(H/\dot{\phi}_{0})\varphi_{\mathbf{k}}$, where
$\dot{\phi}_{0}$ is the rolling speed of the inflaton background.

The approximate scale invariance dictates that the shape function
$\mathcal{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}\rightarrow 0$ and $k_{1}\simeq k_{2}$, it is convenient to choose
the ratio $\varrho=k_{3}/k_{1}$ and the angle $\vartheta$ between
$\mathbf{k}_{1}$ and $\mathbf{k}_{3}$ as the two independent variables, so
that $\mathcal{S}=\mathcal{S}(\varrho,\vartheta)$. The $\vartheta$-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
$\varrho$ 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 $\varrho\rightarrow 0$ limit but can contain
a nonanalytic piece in the form of a noninteger power $\varrho^{\alpha}$ where
$\alpha$ is generally a complex number. Therefore, we expect that the shape
function behaves like (suppressing the angular dependence)

$\boxed{~{}~{}~{}~{}\vspace{3mm}\lim_{\varrho\rightarrow
0}\mathcal{S}(\varrho)=\underbrace{A\varrho^{N}\big{[}1+\mathcal{O}({\varrho})\big{]}}_{\text{analytic}}+\underbrace{B\varrho^{L}\big{[}\sin\big{(}\omega\log\varrho+\varphi\big{)}+\mathcal{O}({\varrho})\big{]}}_{\text{nonanalytic}}.~{}~{}~{}~{}}$
(2)

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 $\omega\neq 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 $\varrho$, described by two
parameters, the frequency $\omega$ and the phase $\varphi$.

It has been noticed that the dimensionless frequency $\omega$ 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
$\omega=\sqrt{m^{2}/H^{2}-9/4}$ when $m>3H/2$. But here we emphasize that all
four parameters ($B$, $L$, $\omega$, $\varphi$) 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\varrho^{N}$ can be thought of as a result of
local (point-like) 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\gg H$). So it corresponds to the
effective field theory (EFT) limit when $m\gg 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\sim$ 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
$\varrho$ 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}\simeq k_{2}\gg k_{3}$, the folded shape $k_{1}\simeq
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}\simeq k_{2}\gg k_{3}$. A typical example of local shape
is

$\mathcal{S}_{\text{local}}\sim\displaystyle\frac{\,{k_{1}^{2}}\,}{k_{2}k_{3}}+\text{2
perms}.$ (3)

This shape is called local because it can arise from a nonlinear field
redefinition
$\zeta(\mathbf{x})\rightarrow\zeta(\mathbf{x})+\lambda\zeta^{2}(\mathbf{x})$
at late universe after inflation, where $\zeta(\mathbf{x})$ is originally a
Gaussian random field and $\lambda$ 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

$\mathcal{S}_{\text{equil}}\sim\displaystyle\frac{\,{k_{1}k_{2}k_{3}}\,}{(k_{1}+k_{2}+k_{3})^{3}}.$
(4)

This shape is often generated by a local interaction in the bulk, such as
$(\varphi^{\prime})^{3}$. (Here a prime denotes conformal time derivative
${\mathrm{d}}/{\mathrm{d}}\tau$; 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 $(\varphi^{\prime})^{3}$ and
$\varphi^{\prime}(\partial_{i}\varphi)^{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.”

Figure 1: The shape function $S$ of the 3-point correlator as a function of
momentum ratio $\varrho=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.

As mentioned above, the nonanalytic piece can be characterized by a
4-parameter set ($B$, $L$, $\omega$, $\varphi$). 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
$\omega$ have simple parametric dependence on the model parameters such as
particle’s mass and spin. The overall amplitude $B$ and the phase $\varphi$,
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 $\mu$) 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,444In most 1-loop processes we have $L=2$. But in the case of a
Dirac fermion running in the loop and without chemical potential enhancement,
the $L=2$ part of the fermion propagator happens to cancel out, and thus the
leading order contribution when $\varrho\rightarrow 0$ appears with $L=3$ [8,
21].

$\displaystyle\lim_{\varrho\rightarrow
0}\mathcal{S}_{\text{signal}}^{\text{(tree)}}(\varrho)\sim\,\text{Re}\,\varrho^{1/2+\mathrm{i}\omega},$
$\displaystyle\lim_{\varrho\rightarrow
0}\mathcal{S}_{\text{signal}}^{\text{(1-loop)}}(\varrho)\sim\,\text{Re}\,\varrho^{2+\mathrm{i}\omega}.$
(5)

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.

Table 1: The parameters of nonanalytic piece in (2) in several representative processes mediated by particles with spin $s$, mass $m$, and chemical potential $\mu$. The value of $B$ only includes leading dependence on $m$ and $\mu$ in the large mass limit. In this table, we set $H=1$.  |  | $B$ |  $L$ | $\omega$
---|---|---|---|---
$s=0$, $m>\frac{3}{2}$, $\mu=0$ [7] | tree | $e^{-\pi m}$ | $\frac{1}{2}$ | $\sqrt{m^{2}-\frac{9}{4}}$
$s=0$, $0<m<\frac{3}{2}$, $\mu=0$ [7] | tree | – | $\frac{1}{2}-\sqrt{\frac{9}{4}-m^{2}}$ | $0$
$s>0$, $m>s-\frac{1}{2}$, $\mu=0$ [13] | tree | $e^{-\pi m}$ | $\frac{1}{2}$ | $\sqrt{m^{2}-(s-\frac{1}{2})^{2}}$
$s>0$, $0<m<s-\frac{1}{2}$, $\mu=0$ [13] | tree | – | $\frac{1}{2}-\sqrt{(s-\frac{1}{2})^{2}-m^{2}}$ | $0$
$s=0$, $m>\frac{3}{2}$, $\mu=0$ [7] | 1-loop | $e^{-2\pi m}$ | $2$ | $2\sqrt{m^{2}-\frac{9}{4}}$
Dirac fermion, $m>0$, $\mu=0$ [8] | 1-loop | $e^{-2\pi m}$ | $3$ | $2m$
Dirac fermion, $m>0$, $\mu>0$ [8] | 1-loop | $e^{2\pi\mu-2\pi\sqrt{m^{2}+\mu^{2}}}$ | $2$ | $2\sqrt{m^{2}+\mu^{2}}$
$s=1$, $m>\frac{1}{2}$, $\mu\geq 0$ [10] | 1-loop | $e^{2\pi\mu-2\pi m}$ | $2$ | $2\sqrt{m^{2}-\frac{1}{4}}$

###  2.2 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 $\zeta$, or equivalently, the
correlators of the inflaton perturbation $\varphi$. Upon quantization, this
correlator can be interpreted as the quantum expectation value of operator
products $\varphi_{\mathbf{k}_{1}}(\tau)\cdots\varphi_{\mathbf{k}_{n}}(\tau)$
at the end of inflation. Throughout the paper, we work with the dS metric
${\mathrm{d}}s^{2}=a^{2}(\tau)(-{\mathrm{d}}\tau^{2}+{\mathrm{d}}\mathbf{x}^{2})$
expressed in the conformal coordinates $(\tau,\mathbf{x})$, and the scale
factor $a(\tau)=1/(-H\tau)$. In particular, the conformal time
$\tau\rightarrow(-\infty,0)$ and thus the end of inflation can be thought of
as the $\tau\rightarrow 0$ limit. We assume the standard Bunch-Davies
condition for the initial state $|\text{BD}\rangle$. So the $n$-point
correlator can be written as
$\langle\varphi_{\mathbf{k}_{1}}\cdots\varphi_{\mathbf{k}_{n}}\rangle$=$\langle\text{BD}|\varphi_{\mathbf{k}_{1}}(\tau_{f})\cdots\varphi_{\mathbf{k}_{n}}(\tau_{f})|\text{BD}\rangle$,
where $\tau_{f}\rightarrow 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):

$\displaystyle\langle\varphi_{\mathbf{k}_{1}}\cdots\varphi_{\mathbf{k}_{n}}\rangle=$
$\displaystyle~{}\int\mathcal{D}\varphi_{+}\mathcal{D}\varphi_{-}\,\varphi_{\mathbf{k}_{1},+}(\tau_{f})\cdots\varphi_{\mathbf{k}_{n},+}(\tau_{f})\,e^{\mathrm{i}S[\varphi_{+}]-\mathrm{i}S[\varphi_{-}]}\delta\big{(}\varphi_{+}(\tau_{f})-\varphi_{-}(\tau_{f})\big{)}.$
(6)

As above, we have two sets of fields $\varphi_{\pm}$, with identical action
$S[\varphi_{\pm}]$, but with an additional minus sign in front of
$S[\varphi_{-}]$ to account for the “backward time”. The two sets of fields
are demanded equal at the future infinity by the $\delta$-function at
$\tau=\tau_{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
$\mathsf{a}=\pm$, corresponding to the two field variables $\varphi_{\pm}$.
The minus type coupling has an additional minus sign in the vertex, coming
from the minus sign in front of $S[\varphi_{-}]$ in the path integral.

Second, each propagator is labeled by two SK indices at the two ends, denoted
as $G_{\mathsf{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 $\varphi_{+}(\tau_{f})=\varphi_{-}(\tau_{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_{\mathsf{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
$\delta$-function, together with an integral over time $\tau$. 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.

##  3 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\pi)^{2}$. In a minimal
scenario, the production is supported by the inflationary expansion and is
efficient for particles with mass $m\lesssim H$. When $m\gg H$, the production
is suppressed by a Boltzmann factor $e^{-\pi 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\pi 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 $\mu\equiv\dot{\phi}_{0}/\Lambda$, where $\Lambda$ is a
cut-off scale and the perturbativity requires $\Lambda>\dot{\phi}_{0}^{1/2}$.
In typical inflation models $\dot{\phi}_{0}^{1/2}\simeq 60H$ so that the
production scale $\mu$ 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 $\mu$. This includes the dim-5 couplings to a
fermion $(\partial_{\mu}\phi)\Psi^{\dagger}\mkern
2.0mu\overline{\mkern-2.0mu\sigma\mkern-2.0mu}\mkern 2.0mu^{\mu}\Psi$ and to a
gauge boson $\phi F\mkern 2.0mu\widetilde{\mkern-2.0muF\mkern-2.0mu}\mkern
2.0mu$. 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
$\mu=\dot{\phi}_{0}/\Lambda$ and this explains why we have a new particle
production mechanism at the scale of $\mu$. 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.

###  3.1 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 $\sigma$ with bare mass
$m_{0}$, which couples to the inflaton through the dim-6 operator
$(\partial_{\mu}\phi)^{2}\sigma^{2}/\Lambda^{2}$. So the Lagrangian is:555 We
adopt the mostly-plus metric throughout the paper.

${\mathscr{L}}=-\sqrt{-g}\bigg{[}\displaystyle\frac{\,{1}\,}{2}(\partial_{\mu}\phi)^{2}+V(\phi)+\displaystyle\frac{\,{1}\,}{2}(\partial_{\mu}\sigma)^{2}+\displaystyle\frac{\,{1}\,}{2}m_{0}^{2}\sigma^{2}-\displaystyle\frac{\,{1}\,}{\Lambda^{2}}(\partial_{\mu}\phi)^{2}\sigma^{2}\bigg{]}.$
(7)

Here $V(\phi)$ is the inflaton potential that generates a slow roll motion
$\dot{\phi}_{0}^{1/2}\simeq 60H\simeq\text{const}$ at the background level.
This rolling generates a mass correction to $\sigma$ so that the $\sigma$ has
an effective mass $m^{2}=m_{0}^{2}+\dot{\phi}_{0}^{2}/\Lambda^{2}$. To avoid
the Boltzmann suppression (see below) we require $m\gtrsim H$. When both terms
in the effective mass are positive, this implies that both terms should be at
most of $\mathcal{O}({H})$. That is,
$\Lambda^{2}\sim\dot{\phi}_{0}^{2}/H^{2}\simeq(60H)^{2}$. The resulting signal
from this choice of parameter is thus very small since the 1-loop process is
suppressed by $(H/\Lambda)^{4}$.

For the current study, we will ignore the issue of this coupling suppression,
insisting that the effective mass $m$ is of $\mathcal{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
$\dot{\phi}_{0}^{2}/\Lambda^{2}$ large but
$m_{0}^{2}-\dot{\phi}_{0}^{2}/\Lambda^{2}$ remains $\mathcal{O}({H})$. 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
$\dot{\phi}_{0}/\Lambda$ and mass differences being $\mathcal{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 $\sigma$ field is at 1-loop order,

(8)

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 functions666Ref. [11]
introduces a scalar chemical potential of a different type which does not work
at the level of mode functions but enhances the signal through particle
injection at interaction vertices.. So an isotropic chemical potential of the
scalar field can only be viewed as a toy example.

###  3.2 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

$\displaystyle{\mathscr{L}}=\sqrt{-g}\bigg{[}-\displaystyle\frac{\,{1}\,}{4}F_{\mu\nu}F^{\mu\nu}-|{\mathrm{D}}_{\mu}\Sigma|^{2}-\lambda|\Sigma|^{4}-\displaystyle\frac{\,{1}\,}{2}(\partial_{\mu}\phi)^{2}-V(\phi)-\displaystyle\frac{\,{1}\,}{\Lambda^{2}}(\partial_{\mu}\phi)^{2}|\Sigma|^{2}\bigg{]}.$
(9)

That is, we consider a scalar QED with $U(1)$ scalar
$\Sigma=(\sigma+\mathrm{i}\pi)/\sqrt{2}$ that couples to the inflaton $\phi$
through the dim-6 operator introduced in the previous subsection. The
difference from the previous case is that the Higgs field $\Sigma$ now has a
nonzero vacuum expectation value
$v^{2}\equiv\langle\sigma\rangle^{2}=\dot{\phi}_{0}^{2}/(\lambda\Lambda^{2})$.
Consequently, there will be a 2-point mixing between $\sigma$ and the inflaton
fluctuation $\varphi$. This is essentially the model considered in [10] except
that we do not include the dim-5 operator $\phi F\mkern
2.0mu\widetilde{\mkern-2.0muF\mkern-2.0mu}\mkern 2.0mu$ here. We will include
this dim-5 operator later. One can also include a bare mass term for $\Sigma$
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_{\xi}$
gauge condition $G=\nabla_{\mu}A^{\mu}-\xi gv\sigma$ where $g$ is the gauge
coupling and $\xi$ 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

$\displaystyle{\mathscr{L}}=$
$\displaystyle~{}\displaystyle\frac{\,{1}\,}{2}A^{\mu}\square
A_{\mu}-\displaystyle\frac{\,{1}\,}{2}a^{2}(g^{2}v^{2}+3H^{2})A^{\mu}A_{\mu}-a^{2}(\partial_{\mu}\bar{c})(\partial^{\mu}c)-a^{4}g^{2}v^{2}\bar{c}c$
$\displaystyle-\displaystyle\frac{\,{1}\,}{2}a^{2}\Big{[}(\partial_{\mu}\sigma)^{2}+(\partial_{\mu}\pi)^{2}+(\partial\varphi)^{2}\Big{]}-\displaystyle\frac{\,{1}\,}{2}a^{4}m_{\sigma}^{2}\sigma^{2}-\displaystyle\frac{\,{1}\,}{2}a^{4}g^{2}v^{2}\pi^{2}$
$\displaystyle-a^{4}\lambda
v\sigma(\sigma^{2}+\pi^{2})-\displaystyle\frac{\,{1}\,}{4}a^{4}\lambda(\sigma^{2}+\pi^{2})^{2}+\displaystyle\frac{\,{\dot{\phi}_{0}v}\,}{\Lambda^{2}}a^{3}\sigma\varphi^{\prime}$
$\displaystyle+a^{2}gA^{\mu}(\sigma\partial_{\mu}\pi-\pi\partial_{\mu}\sigma)-a^{2}g^{2}v\sigma
A^{2}-\displaystyle\frac{\,{1}\,}{2}a^{2}g^{2}(\sigma^{2}+\pi^{2})A^{2}-\displaystyle\frac{\,{1}\,}{2}a^{4}g^{2}v\sigma\bar{c}c+\cdots,$
(10)

where $c$ is the ghost field. We have taken the Feynman-’t Hooft gauge $\xi=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 $\eta_{\mu\nu}$, except in the kinetic term of the gauge
field, where $\square=g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$.

Then at the 1-loop level, we can consider the gauge boson signal from the
following two diagrams.777 Here we only consider diagrams in which the gauge
bosons couple only to the Higgs $\sigma$. There is another class of diagrams
arising from $\phi F\mkern 2.0mu\widetilde{\mkern-2.0muF\mkern-2.0mu}\mkern
2.0mu$ coupling. These diagrams give rise to similar oscillatory signals but
different angular dependence. In [10] it was shown that they are more
suppressed than the ones considered here. For this reason we do not include
these diagrams in the current paper.

(11)

Here we only show the gauge boson diagrams. There are corresponding diagrams
with $\pi$ and ghost fields. We can drop some of these diagrams by considering
special limits. For example, when $\lambda\gg g^{2}$, we only need to include
$\pi$-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 $\lambda\ll g^{2}$ and $v\gg H$. In this limit, we can
drop all graphs containing $\pi$, 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,

$\Delta{\mathscr{L}}=\displaystyle\frac{\,{1}\,}{\Lambda_{F}}\phi F\mkern
2.0mu\widetilde{\mkern-2.0muF\mkern-2.0mu}\mkern 2.0mu.$ (12)

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
$\mu=\dot{\phi}_{0}/\Lambda_{F}$. We refer readers to [10] for more details.

##  4 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 1-loop integrals and
then discuss their numerical implementations. Finally, we present the
numerical results.

###  4.1 One-Loop Integral

Figure 2: The pinched diagram. $\mathsf{a}$ and $\mathsf{b}$ are SK indices.
$\mathbf{k}_{1,2,3}$ are external 3-momenta. $\tau_{1},\tau_{2}$ are conformal
time variables. The black and blue lines represent the boundary-to-bulk
propagators and the internal loop propagators respectively.

All the diagrams from the last subsection, including the triangle diagram
after pinched approximation, can be put into the following form

$\displaystyle\mathcal{B}(k_{1},k_{2},k_{3})=$
$\displaystyle~{}\displaystyle\frac{\,{1}\,}{2}\sum_{\mathsf{a},\mathsf{b}=\pm}\mathsf{ab}(\mathrm{i})^{2}\int\displaystyle\frac{\,{{\mathrm{d}}\tau_{1}}\,}{|H\tau_{1}|^{\alpha}}\displaystyle\frac{\,{{\mathrm{d}}\tau_{2}}\,}{|H\tau_{2}|^{\beta}}\partial_{\tau_{1}}G_{\mathsf{a}}(k_{1};\tau_{1})\partial_{\tau_{1}}G_{\mathsf{a}}(k_{2};\tau_{1})\partial_{\tau_{2}}G_{\mathsf{b}}(k_{3};\tau_{2})\mathcal{I}_{\mathsf{a}\mathsf{b}}^{(1)}(k_{3};\tau_{1},\tau_{2})$
$\displaystyle~{}+(k_{1}\leftrightarrow k_{3})+(k_{2}\leftrightarrow k_{3}),$
(13)

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
$\mathcal{S}(k_{1},k_{2},k_{3})=k_{1}^{2}k_{2}^{2}k_{3}^{2}\mathcal{B}(k_{1},k_{2},k_{3})$.
As explained above, we need to assign SK indices $\mathsf{a}$ and $\mathsf{b}$
to the two vertices, respectively. The factor $|H\tau_{1}|^{\alpha}$ and
$|H\tau_{2}|^{\beta}$ in (13) come from the $\sqrt{-g}$ factor in the
Lagrangian as well as additional $g_{\mu\nu}$ factors in the vertices. We use
$G$ to denote the propagator of the external (inflaton) lines:

$G_{\pm}(k;\tau)=\displaystyle\frac{\,{H^{2}}\,}{2k^{3}}(1\mp\mathrm{i}k\tau)e^{\pm\mathrm{i}k\tau}.$
(14)

In all cases, we include the derivative coupling to the inflaton by using
$\partial_{\tau}G_{\mathsf{a}}(k;\tau)$ for the three external legs. For the
scalar model, there is in principle a contribution from spatial derivative in
the form of $\partial_{i}G\partial^{i}G$ at $\tau_{1}$-vertex. We will drop
this term for simplicity although the contribution from
$\partial_{i}G\partial^{i}G$ is of the same order as
$\partial_{\tau}G\partial_{\tau}G$. The loop integrals
$\mathcal{I}_{\mathsf{ab}}^{(1)}$ are given by

$\mathcal{I}_{\mathsf{a}\mathsf{b}}^{(1)}(k;\tau_{1},\tau_{2})=\int\displaystyle\frac{\,{{\mathrm{d}}^{3}\mathbf{q}}\,}{(2\pi)^{3}}D_{\mathsf{a}\mathsf{b}}(q;\tau_{1},\tau_{2})D_{\mathsf{b}\mathsf{a}}(p;\tau_{2},\tau_{1}).$
(15)

where $D_{\mathsf{a}\mathsf{b}}(k;\tau_{1},\tau_{2})$ are the bulk-to-bulk
propagators

$\displaystyle D_{++}(k;\tau_{1},\tau_{2})=$
$\displaystyle~{}D_{>}(k;\tau_{1},\tau_{2})\theta(\tau_{1}-\tau_{2})+D_{<}(k;\tau_{1},\tau_{2})\theta(\tau_{2}-\tau_{1}),$
(16a) $\displaystyle D_{+-}(k;\tau_{1},\tau_{2})=$
$\displaystyle~{}D_{<}(k;\tau_{1},\tau_{2}),$ (16b) $\displaystyle
D_{-+}(k;\tau_{1},\tau_{2})=$ $\displaystyle~{}D_{>}(k;\tau_{1},\tau_{2}),$
(16c) $\displaystyle D_{--}(k;\tau_{1},\tau_{2})=$
$\displaystyle~{}D_{<}(k;\tau_{1},\tau_{2})\theta(\tau_{1}-\tau_{2})+D_{>}(k;\tau_{1},\tau_{2})\theta(\tau_{2}-\tau_{1}),$
(16d)

with $D_{>}(k;\tau_{1},\tau_{2})\equiv u(\tau_{1},k)u^{*}(\tau_{2},k)$ and
$D_{<}(k;\tau_{1},\tau_{2})\equiv u^{*}(\tau_{1},k)u(\tau_{2},k)$ where
$u(\tau,k)$ is the mode function for the loop particle, either scalar or gauge
boson. It turns out that the four loop integrals ($\mathcal{I}_{++}^{(1)}$,
$\mathcal{I}_{+-}^{(1)}$, $\mathcal{I}_{-+}^{(1)}$, and
$\mathcal{I}_{--}^{(1)}$) are all related to the following unique integral by
including various time ordering and complex conjugation:

$\displaystyle\mathcal{I}_{-+}^{(1)}(k;\tau_{1},\tau_{2})=\int\displaystyle\frac{\,{{\mathrm{d}}^{3}\mathbf{q}}\,}{(2\pi)^{3}}D_{>}(q;\tau_{1},\tau_{2})D_{>}(p;\tau_{1},\tau_{2}).$
(17)

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 $\mu$ to
enhance the signal. In this case, we have $\alpha=2$, $\beta=3$ where
$\alpha,\beta$ are defined in (13), and

$D_{>}(k;\tau_{1},\tau_{2})=\frac{e^{\pi\tilde{\mu}}}{2k}|H\tau_{1}||H\tau_{2}|\text{W}_{-\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(2\mathrm{i}k\tau_{1})\text{W}_{\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(-2\mathrm{i}k\tau_{2}),$
(18)

where W is the Whittaker W function, $\mkern
2.0mu\widetilde{\mkern-2.0mu\mu\mkern-2.0mu}\mkern 2.0mu\equiv\mu/H$ and
$\mkern 2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu\equiv\sqrt{(m/H)^{2}-9/4}$. We only consider scalars with mass $m>3H/2$
so $\mkern 2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern 2.0mu$ is always
real.

Gauge boson loop without chemical potential.

For this example, we will stay in the Feynman-’t Hooft gauge ($\xi=1$). We
only consider the limit $\lambda\ll g^{2}$ and $v\gg 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
$\alpha=0$, $\beta=1$. The propagator will be taken to be,

$\displaystyle D_{\mu\nu>}(k;\tau_{1},\tau_{2})={}$ $\displaystyle
g_{\mu\nu}D_{>}^{\text{(spin-1)}}(k;\tau_{1},\tau_{2})=\frac{g_{\mu\nu}}{2k}\text{W}_{0,\mathrm{i}\tilde{\nu}}(2\mathrm{i}k\tau_{1})\text{W}_{0,\mathrm{i}\tilde{\nu}}(-2\mathrm{i}k\tau_{2})$
$\displaystyle={}$ $\displaystyle\frac{g_{\mu\nu}\pi
e^{-\pi\tilde{\nu}}(\tau_{1}\tau_{2})^{1/2}}{4}\text{H}_{\mathrm{i}\tilde{\nu}}^{(1)}(-k\tau_{1})\text{H}^{(2)}_{-\mathrm{i}\tilde{\nu}}(-k\tau_{2}),$
(19)

where $\text{H}^{(1)}$ ($\text{H}^{(2)}$) is the Hankel function of the first
(second) kind. For the gauge boson of spin-1, we have $\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu=\sqrt{(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;\tau_{1},\tau_{2})D_{>}(p;\tau_{1},\tau_{2})\rightarrow
D_{\mu\nu>}(q;\tau_{1},\tau_{2})D_{>}^{\mu\nu}(p;\tau_{1},\tau_{2})$. Again
the indices here are raised by flat Minkowski metric. The propagator being
proportional to $g_{\mu\nu}$ is a consequence of choosing $\xi=1$ gauge. This
naively introduces $g_{\mu\nu}g^{\mu\nu}=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
$D_{\mu\nu>}(q;\tau_{1},\tau_{2})D_{>}^{\mu\nu}(p;\tau_{1},\tau_{2})\rightarrow
3D^{\text{(spin-1)}}_{>}(q;\tau_{1},\tau_{2})D^{\text{(spin-1)}}_{>}(p;\tau_{1},\tau_{2})$.

Gauge boson loop with chemical potential.

In this case, we still have $\alpha=0$, $\beta=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

$\displaystyle D_{\mu\nu>}(k;\tau_{1},\tau_{2})={}$
$\displaystyle\sum_{h}e^{(h)}_{\mu}(\mathbf{k})e^{(h)*}_{\nu}(\mathbf{k})D_{>}^{(h)}(k;\tau_{1},\tau_{2})\approx
e_{\mu}^{-}(\mathbf{k})e^{-*}_{\nu}(\mathbf{k})D_{>}^{-}(k;\tau_{1},\tau_{2})$
$\displaystyle={}$ $\displaystyle
e_{\mu}^{-}(\mathbf{k})e^{-*}_{\nu}(\mathbf{k})\frac{e^{\pi\tilde{\mu}}}{2k}\text{W}_{-\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(2\mathrm{i}k\tau_{1})\text{W}_{\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(-2\mathrm{i}k\tau_{2}),$
(20)

where $e^{(h)}(\mathbf{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

$\displaystyle\mathcal{I}_{-+}^{(1)}(k;\tau_{1},\tau_{2})=\int\displaystyle\frac{\,{{\mathrm{d}}^{3}\mathbf{q}}\,}{(2\pi)^{3}}D_{\text{reg}}(q;\tau_{1},\tau_{2})D_{\text{reg}}(p;\tau_{1},\tau_{2}),$
(21)

where the regularized propagator $D_{\text{reg}}$ is the original propagator
with mass $m$ subtracted by the same propagator but with a heavier mass $M$,

$D_{\text{reg}}(q;\tau_{1},\tau_{2})\equiv
D_{m}(q;\tau_{1},\tau_{2})-D_{M}(q;\tau_{1},\tau_{2}).$ (22)

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

$\Delta{\mathscr{L}}\supset\delta_{C}(\partial_{\mu}\varphi)^{2}\varphi^{\prime}$
(23)

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 $\delta_{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 $\delta_{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 $\mu$ 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.

###  4.2 Setup and Procedure

We separate the numerical evaluation of the integral (13) into two stages:

  1. (i)

For each combination of SK indices $({\mathsf{a}\mathsf{b}})$, we generate a
discrete 2-dim grid of the conformal time $(\tau_{1},\tau_{2})$ and compute
the loop integral $\mathcal{I}_{\mathsf{a}\mathsf{b}}^{(1)}$ 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 $({\mathsf{a}\mathsf{b}})=(++)$
and $(+-)$. No evaluation for $({\mathsf{a}\mathsf{b}})=(--)$ or $(-+)$ is
needed. Second, as explained in detail later, we only need to evaluate the
$\tau_{1}\geq\tau_{2}$ part of the $({\mathsf{a}\mathsf{b}})=(++)$ and $(+-)$
grids, which reduces about half of the computation task.

  2. (ii)

After obtaining the tabulated grids of
$(\tau_{1},\tau_{2},\mathcal{I}_{\mathsf{a}\mathsf{b}}^{(1)})$ with
$({\mathsf{a}\mathsf{b}})=(++)$ and $(+-)$ from Stage (i), we dress
$\mathcal{I}_{\mathsf{a}\mathsf{b}}^{(1)}$ with the bulk-to-boundary
propagators and other pre-factors as shown in (13). We then interpolate the
dressed grids and integral over $\tau_{1}$ and $\tau_{2}$ to get the 1-loop
result, $\mathcal{B}_{\mathsf{a}\mathsf{b}}(k_{1},k_{2},k_{3})$, for a
particular pair of SK indices $(\mathsf{a}\mathsf{b})$ 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.,
$\mathcal{B}(k_{1},k_{2},k_{3})=\sum_{\mathsf{a},\mathsf{b}=\pm}\mathcal{B}_{\mathsf{a}\mathsf{b}}(k_{1},k_{2},k_{3})+\text{2
perms}=2\text{Re}[\mathcal{B}_{++}(k_{1},k_{2},k_{3})+\mathcal{B}_{+-}(k_{1},k_{2},k_{3})]+\text{2
perms}$, 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., $\text{W}_{\kappa,\lambda}(z)$ with
$\kappa,\lambda\in\mathbb{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,
$\text{W}_{\mp\mathrm{i}\tilde{\mu},\pm\mathrm{i}\tilde{\nu}}(\pm
2\mathrm{i}k\tau)$, are highly oscillatory when $\mkern
2.0mu\widetilde{\mkern-2.0mu\mu\mkern-2.0mu}\mkern 2.0mu$ or $\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern 2.0mu$ is large. As a
result, numerical evaluations for the loop integral
$\mathcal{I}_{\mathsf{a}\mathsf{b}}^{(1)}$, in general, takes a longer time as
$\tilde{\mu}$ or $\tilde{\nu}$ increases. Here we choose $\tilde{\mu}$ and
$\tilde{\nu}$ to be $\mathcal{O}(1)$. For the Pauli-Villars regulator (22), we
set the mass parameter $M$ such that $\tilde{\nu}_{\text{reg}}$ of the
Whittaker functions
$\text{W}_{\mp\mathrm{i}\tilde{\mu},\pm\mathrm{i}\tilde{\nu}_{\text{reg}}}(\pm
2\mathrm{i}k\tau)$ inside $D_{M}(k,\tau_{1},\tau_{2})$, satisfies
$\tilde{\nu}_{\text{reg}}\geq 2\tilde{\nu}$.

  * •

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\equiv
k_{1}=k_{2}\gg k_{3}$), where the CC signal is expected to appear. The other
is a more general range of momenta, the near-equilateral limit
($k_{1}+k_{2}\geq k_{3}$, $k_{1}\leq k_{3}$, $k_{2}\leq k_{3}$), which is
often adopted in literature for reporting the bispectrum. Signals with
$k_{1}=k_{2}\leq k_{3}$ in the near-equilateral limit can be viewed as a
natural extension of those with $k_{1}=k_{2}\gg 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\in[1,k_{\max}]$ with $k_{\max}\sim 10^{3}$, where we expect to see
$\mathcal{O}(\tilde{\nu})$ oscillations per decade in $k$. For the near-
equilateral limit, we consider the triangle region with
$k_{1,2}\in[k_{\min},1]$ with $k_{\min}\sim 0.1$.

  * •

Range and spacing of the time grid. In (13), we take
$\tau_{1,2}\in(-\infty,0)$ for the time integral. But for the numerical
evaluation, the integral range should be finite. Therefore we need to
introduce a finite interval $\tau_{1,2}\in[\tau_{i},\tau_{f}]$. To capture all
relevant physics, the initial time $\tau_{i}$ needs to be set early enough so
that all relevant modes are deep inside the horizon ($|k\tau_{i}|\gg 1$) and
well before the particle production ($|k\tau_{i}|\gg\tilde{\mu}$). The final
time $\tau_{f}$ needs to be late enough such that all modes are well outside
the horizon ($|k\tau_{f}|\ll 1$). Given $k\in[k_{\min},k_{\max}]$ with
$k_{\min}\sim 0.1$ and $k_{\max}\sim 10^{3}$, we choose $\tau_{i}=-200$ and
$\tau_{f}=-10^{-5}$ and log-evenly spacing the interval
$\tau_{1,2}\in[\tau_{i},\tau_{f}]$ into $N_{\text{grid}}=440$ pieces. In the
end, we get a 2-dim time grid of $(\tau_{1},\tau_{2})$ with
$(N_{\text{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|\tau_{i}|-\ln|\tau_{f}|)/N_{\text{grid}}]^{2}$, which we
will use as the measure for the Riemann sum in Stage (ii).

Figure 3: (_Left_) 2-dim time grid in $(\ln|\tau_{1}|,\ln|\tau_{2}|)$ for the
pinched 1-loop diagram with SK indices $(\mathsf{a}\mathsf{b})=(++)$. We only
need to evaluate the loop integrals for the grid with $\tau_{1}\geq\tau_{2}$
(colored part). The loop integral values for the grid with $\tau_{1}<\tau_{2}$
(white part) can be copied from the evaluated grid up to the exchange of
$\tau_{1}\leftrightarrow\tau_{2}$. A similar reduction of computational task
can be achieved for the $(\mathsf{a}\mathsf{b})=(+-)$ diagram, where we only
need to evaluate the grid with $\tau_{1}\geq\tau_{2}$. See text for more
details. (_Right_) An overlay of the integration contour of $\tau$ on top of
the complex plot for the Whittaker function
$\text{W}_{-4\mathrm{i},4\mathrm{i}}(2\mathrm{i}\mathsf{a}\tau)$ with SK index
$\mathsf{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 $\tau\in[\tau_{i},\tau_{f}]$ to
$-\mathrm{i}\tau\in[-\mathrm{i}\tau_{i}-\epsilon,-\mathrm{i}\tau_{f}-\epsilon]$,
where $\epsilon$ is a small positive real number.

  * •

Reducing the grid evaluation for the diagram with
$(\mathsf{a}\mathsf{b})=(++)$. For both scalar and gauge boson loop diagrams,
the key part of the CC signal is given in the form of

$\mathcal{B}_{\mathsf{a}\mathsf{b}}=\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}\int_{-\infty}^{0}{\mathrm{d}}\tau_{2}\int\displaystyle\frac{\,{{\mathrm{d}}^{3}\mathbf{q}}\,}{(2\pi)^{3}}\cdots
D_{\mathsf{a}\mathsf{b}}(q;\tau_{1},\tau_{2})D_{\mathsf{b}\mathsf{a}}(p;\tau_{2},\tau_{1}),$
(24)

where we only keep track terms of interest. For the
$(\mathsf{a}\mathsf{b})=(++)$ diagram, the loop propagator is a piece-wised
function,
$D_{++}(q;\tau_{1},\tau_{2})=D_{>}(q;\tau_{1},\tau_{2})\theta(\tau_{1}-\tau_{2})+D_{<}(q;\tau_{2},\tau_{1})\theta(\tau_{2}-\tau_{1})$,
as shown in Eq. (16). Using the relation
$D_{<}(q;\tau_{1},\tau_{2})=D_{>}(q;\tau_{2},\tau_{1})$, the integral (24) can
be expressed as

$\displaystyle\mathcal{B}_{++}={}$
$\displaystyle\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}\int_{-\infty}^{\tau_{1}}{\mathrm{d}}\tau_{2}\int\displaystyle\frac{\,{{\mathrm{d}}^{3}\mathbf{q}}\,}{(2\pi)^{3}}\cdots
D_{>}(q;\tau_{1},\tau_{2})D_{>}(p;\tau_{1},\tau_{2})$
$\displaystyle+\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}\int^{0}_{\tau_{1}}{\mathrm{d}}\tau_{2}\int\displaystyle\frac{\,{{\mathrm{d}}^{3}\mathbf{q}}\,}{(2\pi)^{3}}\cdots
D_{>}(q;\tau_{2},\tau_{1})D_{>}(p;\tau_{2},\tau_{1})$ (25) $\displaystyle={}$
$\displaystyle\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}\int_{-\infty}^{0}{\mathrm{d}}\tau_{2}\int\displaystyle\frac{\,{{\mathrm{d}}^{3}\mathbf{q}}\,}{(2\pi)^{3}}\cdots$
$\displaystyle\left[D_{>}(q;\tau_{1},\tau_{2})D_{>}(p;\tau_{1},\tau_{2})\theta(\tau_{1}-\tau_{2})+D_{>}(q;\tau_{2},\tau_{1})D_{>}(p;\tau_{2},\tau_{1})\theta(\tau_{2}-\tau_{1})\right],$
(26)

where the two terms inside the square bracket of the last line of (26) are
identical with the exchange of $\tau_{1}\leftrightarrow\tau_{2}$. The same
symmetry can be found when replacing the propagator $D_{++}$ with the
regulated propagator $D_{\text{reg},++}$. Given the symmetry between the two
terms, we only need to evaluate
$\mathcal{I}_{++}^{(1)}(k_{3};\tau_{1},\tau_{2})$ with $\tau_{1}>\tau_{2}$.
After evaluating the loop integrals and obtain the grid
$(\tau_{1},\tau_{2},\mathcal{I}_{++}^{(1)})$ with $\tau_{1}>\tau_{2}$, one can
simply switch $\tau_{1}$ and $\tau_{2}$ and copy the loop integral value
$\mathcal{I}_{++}^{(1)}$ to tabulate the rest of the grid.888Note that
although $\mathcal{I}_{+-}^{(1)}(k_{3},\tau_{1},\tau_{2})$ shares the same
integrand as $\mathcal{I}_{++}^{(1)}(k_{3},\tau_{1},\tau_{2})$ with
$\tau_{1}<\tau_{2}$, we could not simply copy the values between the two given
their differences in the Wick rotation treatment: we rotate
$\tau_{1,2}\rightarrow-\mathrm{i}\tau_{1,2}$
($\tau_{1}\rightarrow+\mathrm{i}\tau_{1}$ and
$\tau_{2}\rightarrow-\mathrm{i}\tau_{2}$) when evaluating
$\mathcal{I}_{++}^{(1)}$ ($\mathcal{I}_{+-}^{(1)}$). See below for Wick
rotation. We illustrate the procedure in the left panel of Fig. 3.

One subtlety arises for the parameter space just around the diagonal axes,
$\tau_{1}=\tau_{2}$. It is the place where the oscillations in the propagators
are significant. If we strictly select grid points that satisfy
$\tau_{1}>\tau_{2}$, the parameter space will be left out given the finite
spacing of the time grid. In the code, we include grid points with
$\tau_{1}=\tau_{2}$, illustrated as darker colored squares in the left panel
of Fig. 3, to represent the values of $\mathcal{I}_{++}^{(1)}$ around them.999
Incidentally, we note that the $\tau_{1}\simeq\tau_{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 $\tau_{1}=\tau_{2}$
region can contribute significantly to the signal as well, depending on the
choice of external momenta.

  * •

Change variables for the time integral. Because the time grids are log-evenly
spaced, it is natural to consider $s_{1,2}\equiv\ln|\tau_{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, $\Delta s_{1}\Delta
s_{2}=[(\ln|\tau_{i}|-\ln|\tau_{f}|)/N_{\text{grid}}]^{2}$. 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.

  * •

Wick rotation of the time variables. The time integral in Stage (ii) contains
a lot of oscillations from the bulk-to-boundary propagators
$\partial_{\tau}G_{\pm}(k,\tau)\sim e^{\pm\mathrm{i}k\tau}$, which is
numerically difficult to handle. Hence we perform the Wick rotation on
$\tau_{1,2}$, for computations in Stage (i) and (ii), to improve numerical
convergence in the early-time limit ($\tau\rightarrow\tau_{i}$).101010There 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. Thanks 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 $\tau\rightarrow+\mathrm{i}\tau$ or
$\tau\rightarrow-\mathrm{i}\tau$, is determined by the
$\mathrm{i}\epsilon$-prescription, which is practically equivalent to such
that the exponential in the bulk-to-boundary propagators are suppressed when
$|\tau|\rightarrow\infty$, i.e., $e^{\pm\mathrm{i}k\tau}\rightarrow
e^{-k|\tau|}=e^{+k\tau}$. Therefore we take $\tau\rightarrow-\mathrm{i}\tau$
($\tau\rightarrow+\mathrm{i}\tau$) if $\tau$ is companied with SK index
$\mathsf{a}=+$ ($\mathsf{a}=-$). In particular, we rotate
$\tau_{1,2}\rightarrow-\mathrm{i}\tau_{1,2}$
($\tau_{1}\rightarrow+\mathrm{i}\tau_{1}$ and
$\tau_{2}\rightarrow-\mathrm{i}\tau_{2}$) when evaluating the loop integral
$\mathcal{I}_{++}^{(1)}$ ($\mathcal{I}_{+-}^{(1)}$).

  * •

Avoiding the branch cut. The Whittaker W function has a branch cut emanating
from 0 to $\infty$. In the implementation of Mathematica 12 and mpmath 1.1,
the branch cut of the Whittaker function $W_{\kappa,\lambda}(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=2\mathrm{i}\mathsf{a}k\tau$ when the SK index $\mathsf{a}=+$. After
performing the Wick rotation described earlier,
$z=2\mathrm{i}\mathsf{a}k\tau=2\mathrm{i}k\tau\rightarrow
z=2k\tau\in\mathbb{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 $2\mathrm{i}k\epsilon$ with
$\epsilon\in\mathbb{R}^{+}$ such that $z=2\mathrm{i}k\tau\rightarrow
z=2k\tau-2\mathrm{i}k\epsilon$ (or equivalently
$\tau\rightarrow-i\tau-\epsilon$). The right panel of Fig. 3 is an
illustration of the rotation procedure for
$\text{W}_{-4\mathrm{i},4\mathrm{i}}(2\mathrm{i}\mathsf{a}\tau)$ with
$\mathsf{a}=+$. Given the range of momentum and time we considered, we fix
$\epsilon=10^{-9}$ and perform the extra shift after rotating $\tau_{1,2}$,
i.e., $\tau_{1,2}\rightarrow-\mathrm{i}\tau_{1,2}-\epsilon$, when evaluating
$\mathcal{I}_{++}^{(1)}$.111111For $\mathcal{I}_{+-}^{(1)}$, the Wick-rotated
$\tau_{1,2}$ do not hit the branch cut under the rotation rule. Hence the
extra small shift is not necessary.

  * •

Reducing the grid evaluation for the diagram with
$(\mathsf{a}\mathsf{b})=(+-)$. Under complex conjugation, the Whittaker
function has the following property

$\left[\text{W}_{-\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(z)\right]^{*}=\text{W}_{+\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(z)\quad\text{for~{}}\tilde{\mu},\tilde{\nu}\in\mathbb{R}~{}\text{and}~{}z\in\mathbb{R}^{+}.$
(27)

We utilize it to reduce the computational load for the
$(\mathsf{a}\mathsf{b})=(+-)$ diagram. The CC signal for the diagram is given
by

$\displaystyle\mathcal{B}_{+-}={}$
$\displaystyle\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}\int_{-\infty}^{0}{\mathrm{d}}\tau_{2}\int\displaystyle\frac{\,{{\mathrm{d}}^{3}\mathbf{q}}\,}{(2\pi)^{3}}\cdots
D_{+-}(q;-\mathrm{i}\tau_{1},\mathrm{i}\tau_{2})D_{-+}(p;\mathrm{i}\tau_{2},-\mathrm{i}\tau_{1})$
$\displaystyle={}$
$\displaystyle\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}\int_{-\infty}^{0}{\mathrm{d}}\tau_{2}\int\displaystyle\frac{\,{{\mathrm{d}}^{3}\mathbf{q}}\,}{(2\pi)^{3}}\cdots
D_{>}(q;\mathrm{i}\tau_{2},-\mathrm{i}\tau_{1})D_{>}(p;\mathrm{i}\tau_{2},-\mathrm{i}\tau_{1}),$
(28)

where we rotated the time variables according to the Wick rotation rules
described above. For each propagator, we found

$\displaystyle
D_{>}(k,\mathrm{i}\tau_{2},-\mathrm{i}\tau_{1};\tilde{\mu},\tilde{\nu})={}$
$\displaystyle\cdots\text{W}_{+\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(-2k\tau_{2})\text{W}_{-\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(-2k\tau_{1})$
$\displaystyle={}$
$\displaystyle\left[\cdots\text{W}_{-\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(-2k\tau_{2})\text{W}_{+\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(-2k\tau_{1})\right]^{*}$
$\displaystyle={}$ $\displaystyle
D_{>}(k,\mathrm{i}\tau_{1},-\mathrm{i}\tau_{2};\tilde{\mu},\tilde{\nu})^{*},$
(29)

where $\cdots$ represents pre-factors that are symmetric under
$\tau_{1}\leftrightarrow\tau_{2}$. We apply Eq. (27) in the second line of Eq.
(29). The same conjugation property is valid for $D_{\text{reg},>}$. Using Eq.
(29), the signal can be split into two “symmetric” parts

$\displaystyle\mathcal{B}_{+-}={}$
$\displaystyle\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}\int_{-\infty}^{0}{\mathrm{d}}\tau_{2}\int\displaystyle\frac{\,{{\mathrm{d}}^{3}\mathbf{q}}\,}{(2\pi)^{3}}\cdots[D_{>}(q;\mathrm{i}\tau_{2},-\mathrm{i}\tau_{1})D_{>}(p;\mathrm{i}\tau_{2},-\mathrm{i}\tau_{1})\theta(\tau_{1}-\tau_{2})$
$\displaystyle+D_{>}(q;\mathrm{i}\tau_{2},-\mathrm{i}\tau_{1})D_{>}(p;\mathrm{i}\tau_{2},-\mathrm{i}\tau_{1})\theta(\tau_{2}-\tau_{1})]$
$\displaystyle={}$
$\displaystyle\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}\int_{-\infty}^{0}{\mathrm{d}}\tau_{2}\int\displaystyle\frac{\,{{\mathrm{d}}^{3}\mathbf{q}}\,}{(2\pi)^{3}}\cdots\left\\{D_{>}(q;\mathrm{i}\tau_{2},-\mathrm{i}\tau_{1})D_{>}(p;\mathrm{i}\tau_{2},-\mathrm{i}\tau_{1})\theta(\tau_{1}-\tau_{2})\right.$
$\displaystyle\left.+\left[D_{>}(q;\mathrm{i}\tau_{1},-\mathrm{i}\tau_{2})D_{>}(p;\mathrm{i}\tau_{1},-\mathrm{i}\tau_{2})\right]^{*}\theta(\tau_{2}-\tau_{1})\right\\},$
(30)

where the two terms inside the curly bracket of the integral (30) are
identical up to the exchange of $\tau_{1}\leftrightarrow\tau_{2}$ and the
complex conjugation. The same symmetry can be found when replacing $D_{>}$
with $D_{\text{reg},>}$. As a result, similar to the case for the
$(\mathsf{a}\mathsf{b})=(++)$ diagram, we only need to evaluate
$\tau_{1}\geq\tau_{2}$ part of the time grid for the loop integral
$\mathcal{I}_{+-}^{(1)}$. After evaluating $\mathcal{I}_{+-}^{(1)}$ for
$\tau_{1}\geq\tau_{2}$, one can simply switch the position of $\tau_{1}$ and
$\tau_{2}$ and conjugate the corresponding $\mathcal{I}_{+-}^{(1)}$ 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
$\text{I}^{(1)}(k_{3};t_{1},t_{2})\equiv(2\pi)^{-3}\int{\mathrm{d}}^{3}\mathbf{q}\cdots
D_{>}(q;t_{1},t_{2})D_{>}(p;t_{1},t_{2})$ with
$t_{1,2}=\mathrm{i}|\tau_{1,2}|-\epsilon$ for $(\mathsf{a}\mathsf{b})=(++)$
and $t_{1}=-\mathrm{i}|\tau_{2}|$ and $t_{2}=\mathrm{i}|\tau_{1}|$ for
$(\mathsf{a}\mathsf{b})=(+-)$. For each grid, we need to evaluate
$(N_{\text{grid}}+1)(N_{\text{grid}}+2)/2$ grid points. Together
$(N_{\text{grid}}+1)(N_{\text{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_{\text{grid}}+1)^{2}$.

  * •

Integral of the loop momentum. For the pinched loop diagram, the momentum
running in the loop propagators are $\mathbf{q}$ and
$\mathbf{p}\equiv\mathbf{q}-\mathbf{k}_{3}$ respectively. We have the freedom
to build the coordinate for $\mathbf{q}$ and we choose a spherical coordinate
with the polar direction ($\theta_{\mathbf{q}}=0$) aligning with the direction
of $\mathbf{k}_{3}$. Under such a coordinate, the magnitude of $\mathbf{p}$ is
given by $p=\sqrt{q^{2}-2qk_{3}\cos\theta_{\mathbf{q}}+k_{3}^{2}}$. And the
integration over $\mathbf{q}$ can be expressed in spherical coordinates
$\int{\mathrm{d}}^{3}\mathbf{q}/(2\pi)^{3}\rightarrow\int_{-1}^{+1}{\mathrm{d}}\cos\theta_{\mathbf{q}}\int_{0}^{+\infty}{\mathrm{d}}q\,q^{2}/(2\pi)^{2}$,
where the azimuth angle of $\mathbf{q}$ has been integrated out.

  * •

Patched Whittaker function. It is quite time-consuming to numerically
integrate the loop integrand over $q\in[0,\infty)$. To speed up the
evaluation, we adopt an approximation for the Whittaker functions, dubbed as
the patched Whittaker functions $\mathscr{W}_{\kappa,\lambda}(z)$: we replace
the Whittaker functions $\text{W}_{\kappa,\lambda}(z)$ with its asymptotic
expansion at $z=0$ and $z=\infty$ for $|z|\ll|\lambda|$ and $|z|\gg|\lambda|$
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|\ll|\lambda|$ ($|z|\gg|\lambda|$) and patched it at $|z|=|\lambda|/5$
($|z|=5|\lambda|$) with the full Whittaker function given our parameter
choices, i.e.,

$\mathscr{W}_{\kappa,\lambda}(z)=\left\\{\begin{matrix}\text{W}^{\text{low}}_{\kappa,\lambda}(z)&|z|<|\lambda|/5\\\
\text{W}_{\kappa,\lambda}(z)&|\lambda|/5\leq|z|<5|\lambda|\\\
\text{W}^{\text{high}}_{\kappa,\lambda}(z)&|z|\geq 5|\lambda|,\\\
\end{matrix}\right.$ (31)

where

$\displaystyle\text{W}^{\text{low}}_{\kappa,\lambda}(z)=\frac{z^{\frac{1}{2}-\lambda}}{24}$
$\displaystyle\left\\{\frac{24\Gamma(2\lambda)}{\Gamma\left(-\kappa+\lambda+\frac{1}{2}\right)}\right.$
$\displaystyle+\frac{z\Gamma(2\lambda-2)}{\Gamma\left(-\kappa+\lambda+\frac{1}{2}\right)}\left[48\kappa(\lambda-1)+3z\left(4\kappa^{2}-2\lambda+1\right)+\frac{\kappa
z^{2}\left(4\kappa^{2}-6\lambda+5\right)}{2\lambda-3}\right]$
$\displaystyle+\frac{z^{2\lambda}\Gamma(-2\lambda-3)}{\Gamma\left(-\kappa-\lambda+\frac{1}{2}\right)}\bigg{[}-48\left(4\lambda^{3}+12\lambda^{2}+11\lambda+3\right)+48\kappa
z\left(2\lambda^{2}+5\lambda+3\right)$
$\displaystyle\left.\left.\phantom{\frac{\Gamma}{\Gamma}}-3z^{2}(2\lambda+3)\left(4\kappa^{2}+2\lambda+1\right)+\kappa
z^{3}\left(4\kappa^{2}+6\lambda+5\right)\right.\bigg{]}\right\\}$ (32)

and

$\text{W}^{\text{high}}_{\kappa,\lambda}(z)=e^{-\frac{z}{2}}z^{\kappa}\left[\frac{-8\lambda^{2}\left(4\kappa^{2}-8\kappa+5\right)+16\lambda^{4}+\left(4\kappa^{2}-8\kappa+3\right)^{2}}{32z^{2}}+\frac{\lambda^{2}-\kappa^{2}+\kappa-\frac{1}{4}}{z}+1\right].$
(33)

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 $(\mathsf{a}\mathsf{b})=(+-)$ diagram.

  * •

Contributions from the momentum permutations. In the above computation, we
focus on the case where the loop integral
$\mathcal{I}^{(1)}_{\mathsf{a}\mathsf{b}}$ is controlled by a fixed external
momentum $k_{3}=1$. To get the contribution from the momentum configuration
$(k_{1}\leftrightarrow k_{3})$ or $(k_{2}\leftrightarrow 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
$\mathcal{S}_{\mathsf{a}\mathsf{b}}(k_{1},k_{2},k_{3})=k_{1}^{2}k_{2}^{2}k_{3}^{2}\mathcal{B}_{\mathsf{a}\mathsf{b}}(k_{1},k_{2},k_{3})$.
The contributions from the two momentum permutations,
$\mathcal{S}_{\mathsf{a}\mathsf{b}}(k_{1}\leftrightarrow k_{3})$ and
$\mathcal{S}_{\mathsf{a}\mathsf{b}}(k_{2}\leftrightarrow k_{3})$, are
respectively given by

$\mathcal{S}_{\mathsf{a}\mathsf{b}}(k_{1}\leftrightarrow
k_{3})=k_{1}^{2}k_{2}^{2}k_{3}^{2}\left(\frac{k_{3}}{k_{1}}\right)^{6}\mathcal{B}_{\mathsf{a}\mathsf{b}}\left(\frac{k_{3}^{2}}{k_{1}},\frac{k_{2}k_{3}}{k_{1}},k_{3}\right)$
(34)

and

$\mathcal{S}_{\mathsf{a}\mathsf{b}}(k_{2}\leftrightarrow
k_{3})=k_{1}^{2}k_{2}^{2}k_{3}^{2}\left(\frac{k_{3}}{k_{2}}\right)^{6}\mathcal{B}_{\mathsf{a}\mathsf{b}}\left(\frac{k_{1}k_{3}}{k_{2}},\frac{k_{3}^{2}}{k_{2}},k_{3}\right),$
(35)

where we rescale $k_{1}$, $k_{2}$, and $k_{3}$ of
$\mathcal{S}_{\mathsf{a}\mathsf{b}}(k_{1}\leftrightarrow k_{3})$
($\mathcal{S}_{\mathsf{a}\mathsf{b}}(k_{2}\leftrightarrow 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
$\mathcal{S}_{\mathsf{a}\mathsf{b}}(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 $\mathcal{O}(10^{5})$ CPU
hours to get the final CC signal for one pinched loop diagram.

###  4.3 Results

Figure 4: (_Left_) The real part of the loop integral $\mathcal{I}_{++}^{(1)}$
(part of the full grid) for a pinched 1-loop diagram of scalar with external
momentum $k_{3}=1$ and parameter
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(3.9,4,8)$. (_Right_) The
real part of the loop integral $\mathcal{I}_{+-}^{(1)}$ (part of the full
grid) with the same parameter setting.

The 1-loop 3-point correlator is the sum of contributions from all SK
diagrams,
$\mathcal{S}(k_{1},k_{2},k_{3})=\sum_{\mathsf{a},\mathsf{b}=\pm}\mathcal{S}_{\mathsf{a}\mathsf{b}}=2\text{Re}\,\mathcal{S}_{++}+2\text{Re}\,\mathcal{S}_{+-}$,
where we applied the reality condition in the last equality. Given the reality
of the dressing functions and pre-factors, the value of
$\text{Re}\,\mathcal{S}_{++}$ ($\text{Re}\,\mathcal{S}_{+-}$) is determined by
the real part of the loop integral $\mathcal{I}_{++}^{(1)}$
($\mathcal{I}_{+-}^{(1)}$). The left and right panels of Fig. 4 respectively
shows
$\text{Re}\,\mathcal{I}_{++}^{(1)}(k_{3}=1;\mathrm{i}|\tau_{1}|-\epsilon,\mathrm{i}|\tau_{2}|-\epsilon)$
and
$\text{Re}\,\mathcal{I}_{+-}^{(1)}(k_{3}=1;\mathrm{i}|\tau_{1}|,-\mathrm{i}|\tau_{2}|)$
as a function of $\ln|\tau_{1}|$ and $\ln|\tau_{2}|$ (parts of the full grids)
for the 1-loop diagram of scalar with parameter
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(3.9,4,8)$. The
integration results contain oscillatory patterns in both
$\text{Re}\,\mathcal{I}_{++}^{(1)}$ and $\text{Re}\,\mathcal{I}_{+-}^{(1)}$.
For the former case, the oscillations are along the axes of
$|\tau_{1}|=|\tau_{2}|$. And for the latter case, the oscillations appear
perpendicular to the axes of $|\tau_{1}|=|\tau_{2}|$. For both cases, the
oscillations become the most significant in the region around
$\ln|\tau_{1}|=\ln|\tau_{2}|\approx 2$.

Figure 5: (_Left_) The rescaled CC signal,
$(k_{1}/k_{3})\mathcal{S}(k_{1},k_{2},k_{3})$, in the squeezed limit,
$k_{1}=k_{2}\gg k_{3}$, for the 1-loop diagram of scalar with
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(3.9,4,8)$ (blue line).
Contributions from the diagrams with SK indices of the same (opposite) sign
are shown as the orange (green) line. The lower sub-panel shows the filtered
CC signal where its non-oscillatory component is attenuated. (_Right_) The CC
signal from diagrams with SK indices of the same sign,
$\mathcal{S}_{++}+\mathcal{S}_{--}=2\text{Re}\,\mathcal{S}_{++}$, rescaled by
$k_{1}/k_{3}$, in the squeezed limit, $k_{1}=k_{2}\gg k_{3}$, for a pinched
1-loop diagram of scalar with parameter
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(5,4,8)$ (solid blue line)
and parameter $(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(5,4,12)$
(dashed orange line). The lower sub-panel shows filtered signals for the two
sets of parameters, where their non-oscillatory components are attenuated. A
set of gray bands, which are evenly separated by
$\Delta\log_{10}(k_{1}/k_{3})=\pi(\log_{10}e)/\tilde{\nu}$, are included for
all the panels as a guide to the oscillatory pattern.

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}\gg 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=\pi\log_{10}e/\tilde{\nu}$. We
insert a set of gray bands, which are evenly separated by
$\Delta\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, $2\text{Re}\,\mathcal{S}_{++}$
($2\text{Re}\,\mathcal{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 $\omega_{c}=0.05\tilde{\nu}/(\pi\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.

Figure 6: (_Upper left_) The rescaled CC signal,
$(k_{1}/k_{3})\mathcal{S}(k_{1},k_{2},k_{3})$, in the squeezed limit,
$k_{1}=k_{2}\gg k_{3}$, for the 1-loop diagram of gauge boson without chemical
potential, where the parameter is set by
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(0,2,4)$ (blue line).
Contributions from the diagrams with SK indices of the same (opposite) sign
are shown as the orange (green) line. The lower sub-panel shows the filtered
CC signal where its non-oscillatory component is attenuated (same for the
other panels). (_Upper right_) The CC signal for the 1-loop gauge boson with
chemical potential, where the parameter is set by
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(3.9,4,8)$. (_Lower left_)
The CC signal for the same diagram with
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(4,4,8)$. (_Lower right_)
The total CC signal for the same diagram with
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(4,4,12)$ (dashed orange)
compared with that of
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(4,4,8)$ (solid blue). A
set of gray bands, which are evenly separated by
$\Delta\log_{10}(k_{1}/k_{3})=\pi(\log_{10}e)/\tilde{\nu}$, are included for
all the panels as a guide to the oscillatory pattern.

Next, we look into the dependence of the signal on the loop regulator,
$\nu_{\text{reg}}$. The right panel of Fig. 5 shows the resulting signal from
the diagrams with the same-sign SK indices,
$\mathcal{S}_{++}+\mathcal{S}_{--}=2\text{Re}\,\mathcal{S}_{++}$, in the
squeezed limit. The solid blue and dashed orange curves are results with
parameter $(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(5,4,8)$ and
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(5,4,12)$ respectively.
The resulting signals from different values of $\nu_{\text{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.

Figure 7: The CC signal $\mathcal{S}(k_{1},k_{2},k_{3})$ in the near-
equilateral limit for the 1-loop diagram of gauge boson with parmeter
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(0,2,4)$ (_upper left_),
$(3.9,4,8)$ (_upper right_), $(4,4,8)$ (_lower left_), and $(4,4,12)$ (_lower
right_). The dashed lines mark the boundaries for the region
$k_{1,2}|\tau_{i}|\gg\max(\tilde{\mu},1)$ when $k_{3}=1$ and $\tau_{i}=-200$.

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
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(0,2,4)$. The oscillatory
signal is relatively tiny since the chemical potential is absent. (Similar
trends are observed for the scalar loop case.) Contributions from the diagrams
with the opposite-sign SK indices become largely suppressed compared to those
from the same-sign SK indices. Once the chemical potential turns on, because
the gauge boson signal is proportional to $\exp(-2\pi h\tilde{\mu})$ where
$h=+$($-$) for its positive (negative) helicity component, the signal from the
negative helicity becomes dominant for $\tilde{\mu}\gtrsim\mathcal{O}(1)$. In
the upper right panel of Fig. 6, we show the resulting signal in the squeezed
limit for the gauge boson with parameter
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(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/opposite-sign SK indices are
comparable. The lower left panel of Fig. 6 shows the results for a slightly
different parameter choice
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(4,4,8)$, where we found
the signal is generally a factor of $\exp(2\pi\Delta\tilde{\mu})\sim 1.8$
larger than that with
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(3.9,4,8)$. We demonstrate
the $\nu_{\text{reg}}$-dependence in the lower right panel, where the solid
blue and dashed orange lines are the resulting total signals with parameter
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(4,4,8)$ and
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{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 $(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})$ described
earlier. We use the dash lines to mark the boundaries of the regions with
$k_{1,2}|\tau_{i}|\gg\max(\tilde{\mu},1)$,121212To be concrete, the boundaries
are set by $k_{1,2}|\tau_{i}|=5\max(\tilde{\mu},1)$.  where the numerical
results are not handicapped by the finite value of $\tau_{i}$ during the time
integral.

###  4.4 Comparison with Analytic Estimates

Figure 8:  (_Upper_) A comparison between the filtered CC signal (blue) and
its analytical estimation (gray). The filtered CC signal is from the 1-loop
diagram of scalar with parameter
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(3.9,4,8)$. Both signals
are rescaled by $k_{1}/k_{3}$. (_Middle_) A similar comparison for the 1-loop
diagram of gauge boson with parameter
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(3.9,4,8)$. (_Lower_) A
similar comparison for the 1-loop diagram of gauge boson with parameter
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(4,4,8)$ (solid blue) or
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(4,4,12)$ (dashed orange).

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 3-point functions. Explicitly, it was
estimated that the signal part of a 1-loop process mediated by a boson with
mass $\mkern 2.0mu\widetilde{\mkern-2.0mum\mkern-2.0mu}\mkern 2.0mu\equiv
m/H\gg 1$ and chemical potential $\mkern
2.0mu\widetilde{\mkern-2.0mu\mu\mkern-2.0mu}\mkern 2.0mu\equiv\mu/H$ takes the
following form as a function of $\varrho\equiv k_{3}/k_{1}$,

$\displaystyle\lim_{\varrho\rightarrow
0}\mathcal{S}_{\text{signal}}^{\text{(est.)}}(\varrho)\simeq\displaystyle\frac{\,{1}\,}{16\pi^{2}}\mkern
2.0mu\widetilde{\mkern-2.0mu\mu\mkern-2.0mu}\mkern 2.0mu^{y}\mkern
2.0mu\widetilde{\mkern-2.0mum\mkern-2.0mu}\mkern 2.0mu^{z}e^{2\pi(\mkern
2.0mu\widetilde{\mkern-2.0mu\mu\mkern-2.0mu}\mkern 2.0mu-\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu)}\varrho^{2}\sin(2\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern 2.0mu\log\varrho+\varphi),$
(36)

where $\mkern 2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu=\sqrt{\mkern 2.0mu\widetilde{\mkern-2.0mum\mkern-2.0mu}\mkern
2.0mu^{2}-9/4}$ for scalar and $\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern 2.0mu=\sqrt{\mkern
2.0mu\widetilde{\mkern-2.0mum\mkern-2.0mu}\mkern 2.0mu^{2}-1/4}$ for spin-1
gauge bosons. The phase $\varphi$ can be calculated and depends on $\mkern
2.0mu\widetilde{\mkern-2.0mu\mu\mkern-2.0mu}\mkern 2.0mu$ and $\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern 2.0mu$ 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\pi(\mkern
2.0mu\widetilde{\mkern-2.0mu\mu\mkern-2.0mu}\mkern 2.0mu-\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern 2.0mu)}$ correctly but
fails to get the correct power dependence $\mkern
2.0mu\widetilde{\mkern-2.0mu\mu\mkern-2.0mu}\mkern 2.0mu^{y}\mkern
2.0mu\widetilde{\mkern-2.0mum\mkern-2.0mu}\mkern 2.0mu^{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 $\varphi$ to match the peaks and valleys of the analytical
and numerical oscillatory signals for the region $k_{1}/k_{3}\gg 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
smaller frequency than analytical estimates
$\tilde{\omega}=2\tilde{\nu}/\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}\gtrsim 20$, as shown
in Fig. 8.

##  5 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. $\phi F\mkern 2.0mu\widetilde{\mkern-2.0muF\mkern-2.0mu}\mkern 2.0mu$
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.

#### Acknowledgment.

We thank Qianhang Ding, Junwu Huang, Edward W. Kolb, and Yi Wang for useful
discussions. We thank Xingang Chen, Soubhik Kumar, Wayne Hu, Junwu Huang,
Hayden Lee, Gustavo Marques-Tavares, and Yi Wang for useful comments on a
draft of this manuscript. This work was completed with resources provided by
the University of Chicago’s Research Computing Center. We also thank Haipeng
An, Lincoln Bryant, Rob Gardner, Pascal Paschos, and Daneng Yang for providing
computational resources and technical supports at the early stage of this
work. LTW and YZ acknowledge the Aspen Center for Physics for hospitality
during the final phase of this study, which was supported by National Science
Foundation grant PHY-1607611 and partially supported by a grant from the
Simons Foundation. LTW is supported by the DOE grant DE-SC0009924. ZZX is
supported by Tsinghua University Initiative Scientific Research Program. YZ is
supported by the Kavli Institute for Cosmological Physics at the University of
Chicago through an endowment from the Kavli Foundation and its founder Fred
Kavli.

##  Appendix A Crosscheck numerical implementations with the tree-level
diagram for complex scalar in an electric field

Ref. [58] investigated the Schwinger effect of charged complex scalars under
an electric field during inflation. The CC signal for the tree-level diagrams
are given by

$\displaystyle\mathcal{B}(k_{1},k_{2},k_{3};\tilde{\mu},\tilde{\nu})=\sum_{\mathsf{a}\mathsf{b}}\mathsf{a}\mathsf{b}(\mathrm{i})^{2}$
$\displaystyle\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}\int_{-\infty}^{0}{\mathrm{d}}\tau_{2}\,\frac{1}{(-H\tau_{1})^{2}}\frac{1}{(-H\tau_{2})^{3}}$
$\displaystyle\partial_{\tau_{1}}G_{\mathsf{a}}(k_{1},\tau_{1})\partial_{\tau_{1}}G_{\mathsf{a}}(k_{2},\tau_{1})\partial_{\tau_{2}}G_{\mathsf{b}}(k_{3},\tau_{2}){D_{\mathsf{b}\mathsf{a}}(k_{3};\tau_{2},\tau_{1};\tilde{\mu},\tilde{\nu})},$
(37)

which contains time integral over $\tau_{1,2}$ and the Whittaker functions in
the bulk-to-bulk propagator
$D_{\mathsf{b}\mathsf{a}}(k_{3};\tau_{2},\tau_{1})$ where $\tilde{\mu}$ and
$\tilde{\nu}$ are the electric field strength parameter and the mass parameter
respectively.131313 Ref. [58] uses $\kappa$ and $\mu$ to represent the
electric field strength parameter and the mass parameter respectively. They
are related to $\tilde{\mu}$ and $\tilde{\nu}$ by
$\kappa=-\mathrm{i}\tilde{\mu}$ and $\mu=\mathrm{i}\tilde{\nu}$. Numerical
evaluations of the diagram take a relatively short time (for codes with the
full Whittaker and the patched Whittaker functions) while some treatments,
such as the Wick rotation of $\tau_{1,2}$ and avoiding the branch cut of the
Whittaker functions, are similar to those for the loop diagram. In addition,
we have an analytical expression for the signal of the
$(\mathsf{a}\mathsf{b})=(+-)$ diagram to compare with the numerical results,

${\mathcal{B}_{+-}(k_{1},k_{2},k_{3};\tilde{\mu},\tilde{\nu})=\frac{\pi^{2}e^{\pi\tilde{\mu}}\left(16\tilde{\nu}^{4}+40\tilde{\nu}^{2}+9\right)}{1024\,\Gamma(1+\mathrm{i}\tilde{\mu})\,k_{1}k_{2}k_{3}^{4}}\sec^{2}(\pi\lambda)\,_{2}\tilde{F}_{1}\left(-\mathrm{i}\tilde{\nu}+\frac{5}{2},\mathrm{i}\tilde{\nu}+\frac{5}{2};-\mathrm{i}\tilde{\mu}+3;\frac{k_{3}-k_{1}-k_{2}}{2k_{3}}\right)},$
(38)

where ${}_{2}\tilde{F}_{1}$ 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, ${\mathcal{S}}(k_{1},k_{2},k_{3})$, in the squeezed
limit $k_{1}=k_{2}\gg k_{3}$, with the electric field strength and the mass
parameter $(\tilde{\mu},\tilde{\nu})=(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.,
$\mathcal{S}=\sum_{\mathsf{a}\mathsf{b}}\mathcal{S}_{\mathsf{a}\mathsf{b}}(\tilde{\mu},\tilde{\nu})+\sum_{\mathsf{a}\mathsf{b}}\mathcal{S}_{\mathsf{a}\mathsf{b}}(-\tilde{\mu},\tilde{\nu})=\sum\mathcal{S}_{\pm\pm}(\tilde{\mu},\tilde{\nu})+\sum\mathcal{S}_{\pm\mp}(\tilde{\mu},\tilde{\nu})+(\tilde{\mu}\rightarrow-\tilde{\mu})=2\text{Re}[\mathcal{S}_{++}(\tilde{\mu},\tilde{\nu})+\mathcal{S}_{+-}(\tilde{\mu},\tilde{\nu})]+(\tilde{\mu}\rightarrow-\tilde{\mu})$,
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
$\sum\mathcal{S}_{\pm\mp}=2\text{Re}\,\mathcal{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\pi\log_{10}e/\tilde{\nu}\approx 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
$\tilde{\mu}$ increases, the overall size of the signal increases, and the
oscillatory patterns of the signals begin to dominate over the smooth
backgrounds.

Figure 9: (_Row 1_) The CC signal $\mathcal{S}$, rescaled by $k_{1}/k_{3}$, in
the squeezed limit, $k_{1}=k_{2}\gg k_{3}$, for a tree-level diagram for
complex scalar in an electric field, where
$\mathcal{S}=2\text{Re}[\mathcal{S}_{++}(\tilde{\mu},\tilde{\nu})+\mathcal{S}_{+-}(\tilde{\mu},\tilde{\nu})+\mathcal{S}_{++}(-\tilde{\mu},\tilde{\nu})+\mathcal{S}_{+-}(-\tilde{\mu},\tilde{\nu})]$.
The electric field strength and the mass parameter are set such that
$(\tilde{\mu},\tilde{\nu})=(\pm 0.1,4),(\pm 1,4),(\pm 4,4)$ for _column 1-3_
respectively. We compare numerical results from the full Whittaker functions
(blue lines) and those from the patched Whittaker functions (orange lines).
(_Row 2-5_) The four components that contribute to $\mathcal{S}$. We also
include the analytical results (dashed black) for $\mathcal{S}_{+-}$
components to compare with the numerical results. A set of gray bands, which
are evenly separated by
$\Delta\log_{10}(k_{1}/k_{3})=2\pi(\log_{10}e)/\tilde{\nu}$, are included for
all panels as a guide to the oscillatory pattern.

###  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
$(\tilde{\mu},\tilde{\nu})=(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 $2\text{Re}\,\mathcal{S}_{+-}$ to the lower
panel of Fig. 10, which agrees well with the numerical results.

Figure 10: (_Upper_) The CC signal $\mathcal{S}$, rescaled by $k_{1}/k_{3}$,
in the squeezed limit, $k_{1}=k_{2}\gg 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 $(\tilde{\mu},\tilde{\nu})=(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_) $\sum\mathcal{S}_{\pm\pm}$ component of $\mathcal{S}$
($\mathcal{S}=\sum\mathcal{S}_{\pm\pm}+\sum\mathcal{S}_{\pm\mp}=2\text{Re}(\mathcal{S}_{++}+\mathcal{S}_{+-})$).
(_Lower_) $\sum\mathcal{S}_{\pm\mp}$ component of $\mathcal{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
$\Delta\log_{10}(k_{1}/k_{3})=2\pi(\log_{10}e)/\tilde{\nu}$, are included for
all panels as a guide to the oscillatory pattern.

##  Appendix B The behavior of the loop integrand at large loop momentum

For the scalar or vector loop integral
$\mathcal{I}^{(1)}_{\mathsf{a}\mathsf{b}}$, its integrand can be expressed in
terms of a single function $\mathcal{L}$ as

$\displaystyle\mathcal{I}^{(1)}(k_{3};\tau_{A},\tau_{B})={}$
$\displaystyle\int_{-1}^{+1}{\mathrm{d}}c_{q}\int_{0}^{\infty}\frac{{\mathrm{d}}q}{(2\pi)^{2}}\mathcal{L}(k_{3};\tau_{A},\tau_{B}),$
(39) $\displaystyle\mathcal{L}(k_{3};\tau_{A},\tau_{B})\equiv{}$
$\displaystyle
q^{2}\tau_{A}^{\alpha}\tau_{B}^{\alpha}\frac{e^{2\tilde{\mu}\pi}}{{4pq}}\text{W}_{-\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(2p\mathrm{i}\tau_{A})\text{W}_{\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(-2p\mathrm{i}\tau_{B})\text{W}_{-\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(2q\mathrm{i}\tau_{A})\text{W}_{\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(-2q\mathrm{i}\tau_{B}),$
(40)

where $c_{q}\equiv\cos\theta_{\mathbf{q}}$,
$p=\sqrt{q^{2}-2qk_{3}c_{q}+k_{3}^{2}}$, and the power $\alpha=2$ ($\alpha=0$)
for the scalar (gauge boson) loop. The loop integral from the
$(\mathsf{a}\mathsf{b})=(++)$ diagram is given by

$\mathcal{I}^{(1)}_{++}(k_{3};\tau_{1},\tau_{2})=\int_{-1}^{+1}{\mathrm{d}}c_{q}\int_{0}^{\infty}\frac{{\mathrm{d}}q}{(2\pi)^{2}}\left[\left.\mathcal{L}(k_{3};\tau_{1},\tau_{2})\right|_{\tau_{1}>\tau_{2}}+\left.\mathcal{L}(k_{3};\tau_{2},\tau_{1})\right|_{\tau_{2}>\tau_{1}}\right],$
(41)

Note that the splitting of the $\tau_{1}>\tau_{2}$ and $\tau_{2}>\tau_{1}$
parts of $\mathcal{L}$ is due to the $\theta$-functions inside the $D_{++}$
propagator. As discussed in Sec. 4.2, we only need to perform the integral for
the first (or the second) $\mathcal{L}$ term inside the integral (41) and the
other part can be tabulated by copying. The loop integral from the
$(\mathsf{a}\mathsf{b})=(+-)$ diagram is given by

$\mathcal{I}^{(1)}_{+-}(k_{3};\tau_{1},\tau_{2})=\int_{-1}^{+1}{\mathrm{d}}c_{q}\int_{0}^{\infty}\frac{{\mathrm{d}}q}{(2\pi)^{2}}\mathcal{L}(k_{3};\tau_{2},\tau_{1}).$
(42)

Again we can split $\mathcal{L}$ into the $\tau_{1}>\tau_{2}$ part and the
$\tau_{1}<\tau_{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 $(\tau_{1},\tau_{2})$. Therefore,
we need to check whether the loop integrand vanishes at the large loop
momentum limit, $q\rightarrow\infty$. Under such limit, the other momentum
inside the loop integral $p\approx q$ since $q-k_{3}\leq p\leq q+k_{3}$. (We
fixed $k_{3}=1$.) Meanwhile, the Whittaker function can be expanded as
$\text{W}_{\kappa,\lambda}(z)\approx
z^{\kappa}e^{-z/2}\left[1+\mathcal{O}(z^{-1})\right]$ at $z\rightarrow\infty$.
Under such limit, Eq. (40) becomes

$\lim_{q\rightarrow\infty}\mathcal{L}(k_{3};\tau_{1},\tau_{2})=\frac{\tau_{1}^{\alpha}\tau_{2}^{\alpha}}{4}\left(\frac{\tau_{1}}{\tau_{2}}\right)^{-2\mathrm{i}\tilde{\mu}}e^{2\mathrm{i}q(\tau_{2}-\tau_{1})}+\mathcal{O}(q^{-1}).$
(43)

Let us first consider the large $q$ behavior of the integrand of
$\mathcal{I}_{++}^{(1)}$. The first $\mathcal{L}$ term inside the integral
(41) is approximately given by

$\displaystyle\lim_{q\rightarrow\infty}\left.\mathcal{L}(k_{3};\tau_{1},\tau_{2})\right|_{\tau_{1}>\tau_{2}}\approx{}$
$\displaystyle\left.\frac{\tau_{1}^{\alpha}\tau_{2}^{\alpha}}{4}\left(\frac{\tau_{1}}{\tau_{2}}\right)^{-2\mathrm{i}\tilde{\mu}}e^{2\mathrm{i}q(\tau_{2}-\tau_{1})}\right|_{\tau_{1}>\tau_{2}}.$
(44)

Because the SK indices associated with $\tau_{1,2}$ are both $+$, we Wick-
rotated $\tau_{1,2}\rightarrow-\mathrm{i}\tau_{1,2}-\epsilon$ where a small
positive number $\epsilon$ is added to avoid the branch cut of the Whittaker
function at $z\in\mathbb{R}^{-}$. Under the Wick-rotation, Eq. (44) yields

$\lim_{q\rightarrow\infty}\left.\mathcal{L}(k_{3};\tau_{1},\tau_{2})\right|_{\tau_{1}>\tau_{2}}\approx\left.\frac{(-\tau_{1}\tau_{2})^{\alpha}}{4}\left(\frac{\tau_{1}}{\tau_{2}}\right)^{-2\mathrm{i}\tilde{\mu}}e^{-2q(|\tau_{2}|-|\tau_{1}|)}\right|_{|\tau_{2}|>|\tau_{1}|}=0$
(45)

given the suppression from $e^{-2q(|\tau_{2}|-|\tau_{1}|)}$. For the same
reasoning the second $\mathcal{L}$ term of the integral (41),
$\left.\mathcal{L}(k_{3};\tau_{2},\tau_{1})\right|_{\tau_{2}>\tau_{1}}$, also
vanishes at large $q$. Next, we consider the large $q$ behavior of the
integrand of $\mathcal{I}_{+-}^{(1)}$. It is given by Eq. (43) up to the
exchange of $\tau_{1}\leftrightarrow\tau_{2}$. Given the SK indices associated
with $\tau_{1}$ and $\tau_{2}$ are $+$ and $-$ respectively, we Wick rotate
$\tau_{1}\rightarrow-\mathrm{i}\tau_{1}$ and
$\tau_{2}\rightarrow+\mathrm{i}\tau_{2}$ and have

$\lim_{q\rightarrow\infty}\mathcal{L}(k_{3};\tau_{2},\tau_{1})\approx\frac{(\tau_{1}\tau_{2})^{\alpha}}{4}\left(-\frac{\tau_{2}}{\tau_{1}}\right)^{-2\mathrm{i}\tilde{\mu}}e^{-2q\left(|\tau_{1}|+|\tau_{2}|\right)}=0$
(46)

given the suppression from $e^{-2q(|\tau_{1}|+|\tau_{2}|)}$.

A tricky point comes from the integrand of $\mathcal{I}_{++}^{(1)}$ at the
$\tau_{1}\rightarrow\tau_{2}$ limit. Under such limit, the exponential
vanishes and no longer provides suppression as $q\rightarrow\infty$. The limit
corresponds to the UV limit of the diagram, i.e., $\Delta\tau\rightarrow
0,q\rightarrow\infty$. As in quantum field theory in flat-space, we replace
the propagators with the regulated propagators, under which the loop integrand
becomes

$\displaystyle\mathcal{L}_{\text{reg}}(k_{3};t_{A},t_{B})={}$ $\displaystyle
q^{2}\tau_{A}^{\alpha}\tau_{B}^{\alpha}\frac{e^{2\tilde{\mu}\pi}}{{4pq}}\left[\text{W}_{-\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(2p\mathrm{i}\tau_{A})\text{W}_{\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(-2p\mathrm{i}\tau_{B})-\text{W}_{-\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}_{\text{reg}}}(2p\mathrm{i}\tau_{A})\text{W}_{\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}_{\text{reg}}}(-2p\mathrm{i}\tau_{B})\right]$
$\displaystyle\left[\text{W}_{-\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(2q\mathrm{i}\tau_{A})\text{W}_{\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}}(-2q\mathrm{i}\tau_{B})-\text{W}_{-\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}_{\text{reg}}}(2q\mathrm{i}\tau_{A})\text{W}_{\mathrm{i}\tilde{\mu},\mathrm{i}\tilde{\nu}_{\text{reg}}}(-2q\mathrm{i}\tau_{B})\right].$
(47)

Expand the Whittaker function to $\text{W}_{\kappa,\lambda}(z)\approx
z^{\kappa}e^{-z/2}\left[1+(\lambda^{2}-\kappa^{2}+\kappa-1/4)z^{-1}+\mathcal{O}(z^{-2})\right]$
under $z\rightarrow\infty$, Eq. (47) at large $q$ is approximated by

$\lim_{q\rightarrow\infty}\mathcal{L}_{\text{reg}}(k_{3};\tau_{1},\tau_{2})\approx\frac{\tau_{1}^{\alpha}\tau_{2}^{\alpha}}{4}\left(\frac{\tau_{1}}{\tau_{2}}\right)^{-2\mathrm{i}\tilde{\mu}}e^{2\mathrm{i}q(\tau_{2}-\tau_{1})}(\tilde{\nu}^{2}_{\text{reg}}-\tilde{\nu}^{2})^{2}\frac{[1-4\tilde{\mu}^{2}+2\tilde{\nu}_{\text{reg}}^{2}+2\tilde{\nu}^{2}+4\mathrm{i}q(\tau_{2}-\tau_{1})]^{2}}{64q^{4}\tau_{1}^{2}\tau_{2}^{2}}.$
(48)

The additional $q^{-4}$ factor provides extra suppressions for the UV limit
(and for integrands with $\tau_{1}\neq\tau_{2}$). Therefore the loop
integrands always vanish at large $q$.

##  Appendix C Crosscheck the filtering methods

In this section, we compare several signal filtering methods that separate the
oscillatory part $O(x\equiv\log_{10}(k_{1}/k_{3}))$ from the non-oscillatory
“background” part $P(x)$ for a signal $S(x)\equiv(k_{1}/k_{3})\mathcal{S}$ in
the squeezed limit. Our first method is the derivative filtering. If $|N|$ and
$|L|$ of (2) is much smaller than $\omega$, $P(x)$ is approximately a
polynomial function of $x$ and vanishes after taking sufficient numbers of
derivatives, while $S(x)$, in the form of $A\sin(\tilde{\omega}x+\varphi)$
with $\tilde{\omega}=\omega/\log_{10}e$, becomes
$A\tilde{\omega}^{n}\sin(\tilde{\omega}x+\varphi^{\prime})$ after taking
$n$-th derivatives. The frequency of the oscillatory pattern $\tilde{\omega}$
keeps the same while its amplitude got enhanced if $\tilde{\omega}>1$. The
derivates unavoidably introduce a phase shift,
$\varphi\rightarrow\varphi^{\prime}$, to the oscillatory pattern. To mitigate
this complication, we take the fourth derivative of $S(x)$ and expect the
oscillatory pattern to be recovered by

$\hat{O}(x)=\tilde{\omega}^{-4}\frac{\partial^{4}}{\partial x^{4}}S(x).$ (49)

The upper panel of Fig. 11 shows an unfiltered signal taken from the lower-
left panel of Fig. 6. It is for a pinched gauge boson 1-loop diagram with
parameter $(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(4,4,8)$ in the
squeezed limit $k_{1}=k_{2}\gg k_{3}$. The signal data is a discrete dataset
with $x$ ranging from $x_{\min}=0$ to $x_{\max}=2$ with an even spacing of
$\Delta x=0.01$. To apply the derivative filtering, we first interpolate the
unfiltered signal. Next, we take the fourth derivatives of the interpolated
function and rescale it with $\tilde{\omega}^{-4}$, where the oscillatory
frequency $\tilde{\omega}=2\tilde{\nu}/\log_{10}e$. The result is shown in the
middle panel of Fig. 11. In practice, one can extract $\tilde{\omega}$ by
fitting $\frac{\partial^{4}}{\partial x^{4}}S(x)$ and checked the pre-
condition $|N|,|L|\ll\omega$ by fitting the resulting $S(x)-\hat{O}(x)$ and
$\hat{O}(x)$.

Figure 11: (_Upper_) The rescaled CC signal
$S=(k_{1}/k_{3})\mathcal{S}(k_{1},k_{2},k_{3})$ in the squeezed limit,
$k_{1}=k_{2}\gg k_{3}$, for the 1-loop diagram of gauge boson with parameter
$(\tilde{\mu},\tilde{\nu},\tilde{\nu}_{\text{reg}})=(4,4,8)$. (_Middle_) The
oscillatory pattern after applying the derivative filtering Eq. (49) to the
signal in the top panel. (_Lower_) The oscillatory patterns after applying
various high-pass filters of Mathematica 12 to the signal in the top panel.
The legend lists the names of the window functions used. A set of gray bands,
which are evenly separated by
$\Delta\log_{10}(k_{1}/k_{3})=\pi(\log_{10}e)/\tilde{\nu}$, are included for
all the panels as a guide to the oscillatory pattern.

An alternative method to access the oscillatory pattern without assuming
particular functional forms for $P(x)$ and $O(x)$ is by using the high-pass
filters. A high-pass filter is a function defined in the frequency space of
the signal. It has a cut-off frequency $\omega_{c}$ and attenuates signals
below the cut-off frequency. The amount of attenuation depends on the filter
window function. Usually, there is a trade-off between the attenuation and so-
called “ringing artifacts”, which refer to the artificial oscillatory signals
introduced by the filtering procedure. A classical example to illustrate this
trade-off is the top-hat filter. The filter is a rectangle function in the
frequency domain, $H(f)=\text{rect}(f)$, and a sinc function in the time
domain, $h(t)=\text{sinc}(t)$. Convoluting the sinc function with a signal
usually results in artificial oscillations. On the opposite end, a Gaussian
filter, which is given by
$G(f)=\exp\left(-\frac{f^{2}}{2\sigma_{f}^{2}}\right)$ in the frequency domain
and
$g(t)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{t^{2}}{2\sigma^{2}}\right)$
with $\sigma=\frac{1}{2\pi\sigma_{f}}$ in the time domain, does not introduce
any ringing artifacts but its power to attenuate low-frequency signal is also
handicapped.

The lower panel of Fig. 11 shows the results after passing the signal through
a set of different high-pass filters that are built-in Mathematica 12
(Hamming, Blackman, Tukey, Hann, and Gaussian) [61]. We set the filter kernel
length to be the length of the signal dataset, and a cut-off frequency
$\omega_{c}=0.147$ for function HighpassFilter for all filter choices. Note
that $\omega_{c}$ needs to satisfy $\omega_{c}\ll\tilde{\omega}/(2\pi)$ to
avoid filtering out $O(x)$ together with $P(x)$. To be concrete, we choose
$\omega_{c}=0.05\,\tilde{\omega}/(2\pi)=0.05\tilde{\nu}/(\pi\log_{10}e)$
through the filtering procedures for the loop processes. The results of
different high-pass filters are broadly consistent with each other. They are
also in reasonable agreement with the result from the derivative filtering.
Therefore we choose the high-pass filter with a Gaussian window function
through the filtering procedures in the main text. Note that for the simple
high-pass filter we adopted, the resulting $O(x)$ may contain a small amount
of the smooth component at small $x$ while the resulting $P(x)=S(x)-O(x)$
contains a small amount of the oscillatory pattern at large $x$. We leave a
detailed study for better filtering procedures for the future.

## References

  * [1] P. D. Meerburg et al., “Primordial Non-Gaussianity,” arXiv:1903.04409 [astro-ph.CO].
  * [2] X. Chen and Y. Wang, “Large non-Gaussianities with Intermediate Shapes from Quasi-Single Field Inflation,” Phys. Rev. D 81 (2010) 063511, arXiv:0909.0496 [astro-ph.CO].
  * [3] X. Chen and Y. Wang, “Quasi-Single Field Inflation and Non-Gaussianities,” JCAP 1004 (2010) 027, arXiv:0911.3380 [hep-th].
  * [4] X. Chen and Y. Wang, “Quasi-Single Field Inflation with Large Mass,” JCAP 1209 (2012) 021, arXiv:1205.0160 [hep-th].
  * [5] S. Pi and M. Sasaki, “Curvature Perturbation Spectrum in Two-field Inflation with a Turning Trajectory,” JCAP 10 (2012) 051, arXiv:1205.0161 [hep-th].
  * [6] J.-O. Gong, S. Pi, and M. Sasaki, “Equilateral non-Gaussianity from heavy fields,” JCAP 11 (2013) 043, arXiv:1306.3691 [hep-th].
  * [7] N. Arkani-Hamed and J. Maldacena, “Cosmological Collider Physics,” arXiv:1503.08043 [hep-th].
  * [8] X. Chen, Y. Wang, and Z.-Z. Xianyu, “Neutrino Signatures in Primordial Non-Gaussianities,” JHEP 09 (2018) 022, arXiv:1805.02656 [hep-ph].
  * [9] L.-T. Wang and Z.-Z. Xianyu, “In Search of Large Signals at the Cosmological Collider,” JHEP 02 (2020) 044, arXiv:1910.12876 [hep-ph].
  * [10] L.-T. Wang and Z.-Z. Xianyu, “Gauge Boson Signals at the Cosmological Collider,” JHEP 11 (2020) 082, arXiv:2004.02887 [hep-ph].
  * [11] A. Bodas, S. Kumar, and R. Sundrum, “The Scalar Chemical Potential in Cosmological Collider Physics,” JHEP 02 (2021) 079, arXiv:2010.04727 [hep-ph].
  * [12] X. Chen, Y. Wang, and Z.-Z. Xianyu, “Loop Corrections to Standard Model Fields in Inflation,” JHEP 08 (2016) 051, arXiv:1604.07841 [hep-th].
  * [13] H. Lee, D. Baumann, and G. L. Pimentel, “Non-Gaussianity as a Particle Detector,” JHEP 12 (2016) 040, arXiv:1607.03735 [hep-th].
  * [14] X. Chen, Y. Wang, and Z.-Z. Xianyu, “Standard Model Background of the Cosmological Collider,” Phys. Rev. Lett. 118 no. 26, (2017) 261302, arXiv:1610.06597 [hep-th].
  * [15] X. Chen, Y. Wang, and Z.-Z. Xianyu, “Standard Model Mass Spectrum in Inflationary Universe,” JHEP 04 (2017) 058, arXiv:1612.08122 [hep-th].
  * [16] H. An, M. McAneny, A. K. Ridgway, and M. B. Wise, “Quasi Single Field Inflation in the non-perturbative regime,” JHEP 06 (2018) 105, arXiv:1706.09971 [hep-ph].
  * [17] S. Kumar and R. Sundrum, “Heavy-Lifting of Gauge Theories By Cosmic Inflation,” JHEP 05 (2018) 011, arXiv:1711.03988 [hep-ph].
  * [18] X. Chen, Y. Wang, and Z.-Z. Xianyu, “Schwinger-Keldysh Diagrammatics for Primordial Perturbations,” JCAP 1712 no. 12, (2017) 006, arXiv:1703.10166 [hep-th].
  * [19] Y.-P. Wu, “Higgs as heavy-lifted physics during inflation,” JHEP 04 (2019) 125, arXiv:1812.10654 [hep-ph].
  * [20] L. Li, T. Nakama, C. M. Sou, Y. Wang, and S. Zhou, “Gravitational Production of Superheavy Dark Matter and Associated Cosmological Signatures,” JHEP 07 (2019) 067, arXiv:1903.08842 [astro-ph.CO].
  * [21] S. Lu, Y. Wang, and Z.-Z. Xianyu, “A Cosmological Higgs Collider,” arXiv:1907.07390 [hep-th].
  * [22] T. Liu, X. Tong, Y. Wang, and Z.-Z. Xianyu, “Probing P and CP Violations on the Cosmological Collider,” arXiv:1909.01819 [hep-ph].
  * [23] A. Hook, J. Huang, and D. Racco, “Searches for other vacua II: A new Higgstory at the cosmological collider,” arXiv:1907.10624 [hep-ph].
  * [24] A. Hook, J. Huang, and D. Racco, “Minimal signatures of the Standard Model in non-Gaussianities,” arXiv:1908.00019 [hep-ph].
  * [25] S. Kumar and R. Sundrum, “Cosmological Collider Physics and the Curvaton,” arXiv:1908.11378 [hep-ph].
  * [26] S. Alexander, S. J. Gates, L. Jenks, K. Koutrolikos, and E. McDonough, “Higher Spin Supersymmetry at the Cosmological Collider: Sculpting SUSY Rilles in the CMB,” JHEP 10 (2019) 156, arXiv:1907.05829 [hep-th].
  * [27] D.-G. Wang, “On the inflationary massive field with a curved field manifold,” JCAP 01 (2020) 046, arXiv:1911.04459 [astro-ph.CO].
  * [28] Y. Wang and Y. Zhu, “Cosmological Collider Signatures of Massive Vectors from Non-Gaussian Gravitational Waves,” arXiv:2001.03879 [astro-ph.CO].
  * [29] L. Li, S. Lu, Y. Wang, and S. Zhou, “Cosmological Signatures of Superheavy Dark Matter,” arXiv:2002.01131 [hep-ph].
  * [30] J. Fan and Z.-Z. Xianyu, “A Cosmic Microscope for the Preheating Era,” JHEP 01 (2021) 021, arXiv:2005.12278 [hep-ph].
  * [31] S. Aoki and M. Yamaguchi, “Disentangling mass spectra of multiple fields in cosmological collider,” JHEP 04 (2021) 127, arXiv:2012.13667 [hep-th].
  * [32] N. Maru and A. Okawa, “Non-Gaussianity from $X,Y$ gauge bosons in Cosmological Collider Physics,” arXiv:2101.10634 [hep-ph].
  * [33] P. D. Meerburg, M. Münchmeyer, J. B. Muñoz, and X. Chen, “Prospects for Cosmological Collider Physics,” JCAP 1703 no. 03, (2017) 050, arXiv:1610.06559 [astro-ph.CO].
  * [34] A. Moradinezhad Dizgah, H. Lee, J. B. Muñoz, and C. Dvorkin, “Galaxy Bispectrum from Massive Spinning Particles,” JCAP 1805 no. 05, (2018) 013, arXiv:1801.07265 [astro-ph.CO].
  * [35] K. Kogai, K. Akitsu, F. Schmidt, and Y. Urakawa, “Galaxy imaging surveys as spin-sensitive detector for cosmological colliders,” JCAP 03 (2021) 060, arXiv:2009.05517 [astro-ph.CO].
  * [36] N. Barnaby and M. Peloso, “Large Nongaussianity in Axion Inflation,” Phys. Rev. Lett. 106 (2011) 181301, arXiv:1011.1500 [hep-ph].
  * [37] N. Barnaby, R. Namba, and M. Peloso, “Phenomenology of a Pseudo-Scalar Inflaton: Naturally Large Nongaussianity,” JCAP 1104 (2011) 009, arXiv:1102.4333 [astro-ph.CO].
  * [38] P. Adshead, E. Martinec, and M. Wyman, “Gauge fields and inflation: Chiral gravitational waves, fluctuations, and the Lyth bound,” Phys. Rev. D88 no. 2, (2013) 021302, arXiv:1301.2598 [hep-th].
  * [39] P. Adshead and E. I. Sfakianakis, “Fermion production during and after axion inflation,” JCAP 11 (2015) 021, arXiv:1508.00891 [hep-ph].
  * [40] M. Peloso, L. Sorbo, and C. Unal, “Rolling axions during inflation: perturbativity and signatures,” JCAP 1609 no. 09, (2016) 001, arXiv:1606.00459 [astro-ph.CO].
  * [41] S. Weinberg, “Quantum contributions to cosmological correlations,” Phys. Rev. D 72 (2005) 043514, arXiv:hep-th/0506236.
  * [42] N. Arkani-Hamed, D. Baumann, H. Lee, and G. L. Pimentel, “The Cosmological Bootstrap: Inflationary Correlators from Symmetries and Singularities,” JHEP 04 (2020) 105, arXiv:1811.00024 [hep-th].
  * [43] D. Baumann, C. Duaso Pueyo, A. Joyce, H. Lee, and G. L. Pimentel, “The cosmological bootstrap: weight-shifting operators and scalar seeds,” JHEP 12 (2020) 204, arXiv:1910.14051 [hep-th].
  * [44] D. Baumann, C. Duaso Pueyo, A. Joyce, H. Lee, and G. L. Pimentel, “The Cosmological Bootstrap: Spinning Correlators from Symmetries and Factorization,” arXiv:2005.04234 [hep-th].
  * [45] C. Sleight, “A Mellin Space Approach to Cosmological Correlators,” JHEP 01 (2020) 090, arXiv:1906.12302 [hep-th].
  * [46] C. Sleight and M. Taronna, “Bootstrapping Inflationary Correlators in Mellin Space,” JHEP 02 (2020) 098, arXiv:1907.01143 [hep-th].
  * [47] E. Pajer, D. Stefanyszyn, and J. Supeł, “The Boostless Bootstrap: Amplitudes without Lorentz boosts,” JHEP 12 (2020) 198, arXiv:2007.00027 [hep-th].