Gauge Boson Signals at the Cosmological Collider

We study the production of massive gauge bosons during inflation from the axion-type coupling to the inflaton and the corresponding oscillatory features in the primordial non-Gaussianity. In a window in which both the gauge boson mass and the chemical potential are large, the signal is potentially reachable by near-future large scale structure probes. This scenario covers a new region in oscillation frequency which is not populated by previously known cosmological collider models. We also demonstrate how to properly include the exponential factor and discuss the subtleties in obtaining power dependence of the gauge boson mass in the signal estimate.


Introduction
Cosmological inflation in the early universe sets the stage for rich dynamics of particle physics at energy scales much above the reach of terrestrial experiments. In the coming decades, much more observational data will further shed light in this era. In particular, the precision in the primordial Non-Gaussianity (NG) measurement will be improved by orders of magnitudes [1]. Among various NG observables, the oscillatory shape in the squeezed limit due to particle production during the inflation is particularly striking. (We will henceforth refer to this oscillatory shape the "signal.") Detecting such a signal at this so-called cosmological collider offers direct evidence of new physics particles and a tool of studying their properties [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21].
The strength of the signal depend sensitively on the coupling between the inflaton and the new physics particles. One key difference between the cosmological collider and a terrestrial collider experiment is that the interaction with the inflaton can change the spectrum of the new physics particles significantly. Very often, the signal size is exponentially sensitive to the mass. Hence, we will only have observable signals with couplings of specific types [19]. This leads us to focus on a specific class of couplings by assuming the inflaton has an approximate shift symmetry, φ → φ + c, which is well motivated by the requirement of slow roll inflation. In [19], it is further argued that a sub-class of such couplings are particularly promising. They are  (47) and (48)]. In this plot we take u = 1 which is defined in (3). Also shown in the figure are the constraints from the validity of EFT expansion (upper grey region), no large back reaction to inflation dynamics (meshed regions), the validity of perturbative calculation (solid and dashed blue curves), and the associated equilateral NG (magenta curve). See Sec. 4 for more discussions.
where J µ5 is a chiral fermion current and F is the field strength of a gauge field. The cut off Λ parameterizes the scale of the physics which is responsible for the generation of these operators. These couplings introduce new sources of particle production during inflation and thus additional enhancement to the signal. In the case of the axial coupling to a fermion of mass m, the chemical potential µ =φ 0 /Λ introduces both an enhancement e πµ/H (through the particle production) and the Boltzmann suppression e −π √ m 2 +µ 2 /H (through the mass correction). Therefore, when µ m, the enhancement can help to largely cancel the Boltzmann suppression, leaving an O(1) rate for particle production. Using this mechanism to generate large signals has been explored in several directions in the context of cosmological collider physics [12,16,17].
The case of gauge boson production is more subtle. In the case where the gauge boson has mass of O(H), the coupling φF F /Λ does not work by itself to generate the signal. On the one hand, a low cutoff Λ would lead to exponentially fast and thus potentially dangerous particle production, which scales as power of e π(µ−m A )/H . On the other hand, to avoid the exponential factor e π(µ−m A )/H , we would want µ =φ/Λ ∼ H which implies a relatively high cutoff Λ ∼φ 0 /H ∼ 3 × 10 3 H. This will suppress the signal size which scales as (H/Λ) 3 . More careful estimate shows that there is no viable intermediate range by varying the cutoff scale alone.
In this paper, we would like to emphasize that the scenario where µ ∼ m A H is much more interesting. First of all, this is a plausible scenario. We in general expect that the heavy gauge boson gets its mass m A = gσ 0 from the background value of some Higgs field σ 0 . It is more natural to consider that σ 0 is linked to dynamics of the inflation. In this case, we expect it to be  [19], and rough reaches of current and future observations. The light gray contours are the predictions of the quasi-single-field inflation scenario and several other scenarios, and the light gray shaded region on the right is the signal of the models in which the inflaton has a coupling to the chiral current of massive fermions. Note that the gauge boson signals could show up in a parameter region (4 ν 20, f (osc) NL 0.1) not populated by previously known scenarios.
higher than H, which is parametrically lower than the inflation scale. For the size of the signal, µ H allows us to have a lower cutoff scale Λ, as long as the EFT bound Λ φ 1/2 0 is satisfied. This will help to relieve the coupling suppression without introducing the large exponential factor. Moreover, the Higgs can mix with the inflaton. Therefore, we have additional graphs generating the signal through the Higgs-gauge coupling and the Higgs-inflaton mixing, which is in general much larger than the graph with the axion coupling alone.
Studying the above signals more carefully is the main purpose of this paper. The main results of this paper are summarized in Figs. 1 and 2. In Fig. 1 we show the signal size of the oscillatory features as functions of gauge boson mass m A and the ratio ≡ µ/m A , together with several constraints. The main conclusion is that there are parameter regions with sizable signals and consistent with current constraints. In Fig. 2 we contrast this signal with signals from previously studied models, together with the expected observation reach. The main message here is that the gauge boson signals studied here could populate an intermediate mass range 5 m/H 20 with sizable signal size, which is in general not possible for other channels.
In [19], we developed a simple way of estimating the sizes of the signal based largely on power counting of the relevant couplings, loop factors, and propagators. For the purpose of studying the signals of gauge boson production, we refined such estimates to include appropriate exponential factors in this paper. We also performed a more detailed study of the behavior of the gauge boson propagator. The main takeaways are 1) It is easy to estimate the parameter dependence correctly for the "EFT part" of the NG, and also easy to estimate correctly the exponential factor of the "signal." However, it is in general hard to get the power dependence in the "signal part." 2) The oscillatory signals are contributed by both the "non-local" and "local" part of the propagator in the late time limit, which was overlooked by many previous studies.
Gauge boson productions through an axionic coupling have been studied in different contexts [22][23][24][25][26][27][28]. In contrast to our paper, most of these studies focused on the scenario in which the gauge boson is massless. Ref. [29] considered the possibility that the gauge symmetry is Higgsed, to achieve better consistency with the CMB power spectrum and gravity wave measurements. In this paper, we will focus on the production of an Abelian gauge boson. There could also be production of non-Abelian gauge bosons which could has additional interesting consequences for the cosmological collider physics. The presence of gauge field could also source a large tensor mode in the primordial fluctuation [30][31][32].
The rest of this paper is organized as follows. In Sec. 2 we provide model motivations for the gauge boson signal. We figure out the signal in Sec. 3 and various constraints in Sec. 4. We conclude in Sec. 5. In App. A we collect some discussions about signal estimate and also about the late-time expansion of the propagators.

Framework
In this section, we layout the framework for our analysis. The starting point is an inflaton endowed with an approximate shift symmetry φ → φ + c. An extensively studied class of inflation models with such a shift symmetry is the axion inflation, which has a long history with many possible scenarios [31,[33][34][35][36][37][38] (see [39] for a review). Motivated by this, we consider the scenario with the following Lagrangian, Here φ is the inflaton and Σ is a complex scalar charged under a local U (1), D µ Σ = (∂ µ + igA µ )Σ. The rolling of the inflaton field ∂ µ φ =φ 0 δ µ0 then generates a vacuum expectation value (VEV) for Σ. Using the parameterization Σ = (σ + iπ)/ √ 2, we can write σ 2 = σ 2 0 = (φ 2 0 /Λ 2 Σ + m 2 Σ )/λ. Here m 2 Σ can have either sign, but we assume that the combinationφ 2 0 /Λ 2 Σ + m 2 Σ is positive so that σ picks up nonzero VEV. Then the scalar mass in this minimum is m 2 σ = 2λσ 2 0 , while the gauge boson mass is m A = gσ 0 . In addition, in the rolling inflaton background, we have a chemical potential to the gauge boson, µ =φ 0 /Λ F .
In general, m A , m σ , and µ are free parameters in this scenario. We will be focusing on the case µ ∼ m A in this paper. To represent the relevant parameter space, we begin with an special limit in which m 2 Σ φ 2 0 /Λ 2 Σ . Hence, we have the VEV σ 2 0 =φ 2 0 /(λΛ 2 Σ ), and m A = gσ 0 = gφ 0 /( √ λΛ Σ ). If we consider g ∼ λ ∼ O(1), both masses m A and m σ will be much higher than the Hubble if the cutoff scale Λ Σ is close to its lower boundφ 1/2 0 . If we further assume Λ Σ = Λ F , the chemical potential µ =φ 0 /Λ F can be close to m A . Of course, we will consider more general cases beyond this simple limit. To this end, we introduce a new parameter u − 1 measures the deviation from the limit m Σ = 0. In general, we could have any u > 0. However, we expect u 1 to be tuned. In general, the cut offs Λ Σ and Λ F could be different. We take this into account by trade Λ F with the chemical potential and treat it as a free parameter. In practice, it is more convenient to define the ratio ≡ µ/m A and use it instead of µ.
Next, we discuss possible models relevant for the signal we study in this paper. Our discussion here is not aiming at constructing a specific model. Instead, it is to motivate the corresponding parameter space of interest through some general consideration and examples. In addition to those mentioned here, there could certainly be other scenarios which can give rise to similar signals.
We begin with the energy scales involved in the problem. The scale of the inflation, Λ inf ≡ ρ 1/4 inf , is constrained by the current observation to be at most ∼ 10 16 GeV. For the following discussion, we will take this upper bound as a benchmark value. This in turn sets the Hubble scale to be The chemical potential is set by a dimension-5 operator in the inflation background, µ =φ 0 /Λ F . From the validity of the EFT expansion, we have Λ F ≥φ 1/2 0 . Hence, µ <φ 1/2 0 60H. For reasons we will discuss in detail, we would focus on the case in which µ ∼ m A . Therefore, we would be mostly interested in considering µ ∼ m A ∼ 10 14∼15 GeV.
Since the inflaton has a shift symmetry, it is natural to consider it as a pseudo-Goldstone boson resulting from some spontaneous symmetry breaking, with the scale f which is also called the decay constant of the inflaton. At the same time, as in many familiar examples, such a Goldstone could be non-linearly realizing a symmetry which has anomalies, resulting in a coupling of the form f eff is an effective decay constant characterizing the specific coupling between the inlfaton φ and the gauge field F . This f eff is to be distinguished from the scale of the symmetry breaking, f . While there is a general expectation f ≤ M Pl [40][41][42][43][44], the decay constant of a generic Goldstone boson can be a free parameter otherwise. However, specializing to the case of the inflaton, it is more natural to consider that f is related to the dynamics around the scale of inflation, and hence not too far away from Λ inf (possibly deferring by small couplings and loop factors). For example, the models from string compactifications considered in [45] have a range f /M Pl > 10 −4 . Hence, a typical benchmark would have f ∼ 10 16 -10 17 GeV. From Eq. (2) and Eq. (4), we have µ ∼φ 0 /(16π 2 f eff ). If we naively take f = f eff , we will have µ ∼ 10 −1 H. Hence, to be in the parameter region relevant to this paper, we would need f ∼ 10 2∼3 f eff . Recently, this has been demonstrated to be achievable in a broad range of models [46,47].
A potential for the inflaton must be generated, making it a pseudo-Goldstone. We would like to check that a successful inflation can happen in the region of parameter space we consider in this paper. Perhaps the easiest way to satisfy such a requirement is to imagine that the decay constant f dose not play a direct role in the inflaton potential. This is the case, for example, for the so-called monodromy motivated models [34,35], which has a potential of the form where W (φ), generated by stringy dynamics, is the dominant piece in the potential which drives the inflation. For example, we could have W (φ) = µ 3 φ, with φ ∼ 10M Pl during inflation and µ ∼ 6 × 10 −4 M Pl . The second term in the potential depends on f and allows the possibility that the axion/inflaton couples to a confining sector. It depends on f . However, with the assumption Λ 4 W , it does not play a significant role in driving inflation. We could also consider more "economical" cases. One such example is the so called Natural Inflation scenario [33], with a potential of the form V (φ) = Λ 4 cos(1 − cos(φ/f inf )), where f inf is an effective decay constant and characterizes the period of the inflaton potential. f inf depends on both the symmetry breaking scale f and the mechanism of generating the potential. However, for the simplest case f inf f < M Pl , this potential is not consistent with observations as it would produce a very red spectrum. A number of mechanisms, typically involving multiple axion-like particles, have been invented to generate an effective decay constant f inf > M Pl > f [36,38,48].
If we still keep f ∼ 10 16 GeV as our benchmark, it is not difficult to imagine that the period of the inflaton can be larger by a factor of 10 3∼4 through one of these mechanisms, and hence this can be a viable inflation model up to potential constraint from weak-gravity arguments [42][43][44]. We will then also invoke one of the mechanisms to enhance the coupling between the inflaton and the spectator gauge field F , as discussed above.
As a concrete example, we could have a scenario in which the inflaton can have the following couplings where a, b ∼ O(1). M, N 1 are numbers which depends on the details of the physics (at a scale ∼ f ) which generates such couplings. For example, in the models presented in Ref. [47], we have N ∼ M 2 . G is the field strength of a gauge field which would become strongly coupled at a scale around Λ inf , and generates a potential for the inflaton φ. In this case, H ∼ Λ 2 inf /M Pl ∼ 10 13 GeV. Hence,φ 1/2 0 60H ∼ 10 15 GeV. For a benchmark value of scale f , we can take f ∼ (several) × 10 16 − 10 17 GeV (> Λ). Therefore, M ∼ 100 gives f inf = M f > M Pl . At the same time, the effective coupling to a spectator gauge field F is f eff = f /N ∼ 10 13 GeV. In this case, the effective chemical potential for gauge field F is µ ∼φ 0 /(16π 2 f eff ) ∼ 10 15 GeV. These are of course rough and parametric estimates, and it is easy to get one order of magnitude either way.
Next, we discuss the mass of the gauge boson of the spectator gauge group. It is easy to imagine it to be Higgsed. In principle, depending on the Higgs potential, its mass can be a free parameter. For the scenario under consideration in the paper, it would be more interesting to consider the case in which the physics at the inflation scale Λ inf is also responsible for the generation of the gauge boson mass. For example, some strong dynamics can generate the inflation scale, such as in the models mentioned above. At the same time, the strong dynamics can also have a richer structure. It can break a global symmetry which give rise to a coset worth of Goldstones, similar to the case of QCD. With additional explicit breaking, these Goldstones can develop a potential which in the end Higgses the spectator gauge group. Such a set up has been explored extensively for the electroweak symmetry breaking, known as the composite Higgs models (see [49] for a review). Implementing a similar mechanism to our setup with the strong coupling scale being Λ inf , we expect the gauge boson mass to be An operator at dimension-6, will also play an important role in our study, both giving an important contribution to the size of the Higgs VEV σ 0 and providing an important coupling to mediate the signal. The size of Λ Σ is determined by the physics which mediates the interaction between the Higgs and the inflaton. Similar to the discussion of other mass scales above, it would be more natural to consider that Λ Σ is also related to the physics which governs the inflation. For example, we can consider the strong coupling example again, where the strong coupling scale is Λ inf ∼ 10 16 GeV. The coupling between a pseudo-Goldstone (Σ) and the inflaton (not part of the composite states) can be mediated by one of the composite resonances with mass m * and coupling g * . In general, we 3 The Signal

Gauge Boson Production
The gauge boson production from the axion coupling has been studied a lot in various contexts. Here we briefly summarize the result relevant to our investigation of cosmological collider signals.
We begin with evaluating the Lagrangian (2) with the scalar background ∂ µ φ =φ 0 δ µ0 and Σ = σ 0 / √ 2, and keeping gauge boson terms only, where again µ =φ 0 /Λ F , m A = gσ 0 . The physical time t and the conformal time τ during inflation are related by e Ht = a = −1/(Hτ ). Imposing the condition ∂ µ ( √ −gA µ ) = 0 removes the A 0 component and then yields the equation of motion for the spatial component A i , which can be written in terms of k-mode as where the prime denotes the conformal time derivative. These three equations can be decoupled by going to the helicity basis Then the decoupled equations read, As explained in [19], the combination µh is a general feature of the chemical potentials that can enhance particle productions during inflation. The dispersion relation for A (±1) is where k phys = k/a is the physical momentum. Working in the limit m A ∼ µ H, we can ignore the last term. From (12), we see that the new feature for the gauge boson's chemical potential, as opposed to the fermionic case, is that the mode equation (11) can become truly tachyonic when µ > m A . As a result, both the positive and the negative frequency part of the mode function in the late-time limit will be exponentially enhanced.
Assuming a Bunch-Davis-like initial condition for A (h) , the solutions to equations (11) are We only consider m A /H 1 in this work so ν is always real and positive. One can readily check that the early-time limit (τ → −∞) has the following properly normalized form, On the other hand, the particle production can be most easily seen by looking at the late-time limit |kτ | 1 of the wavefunction of the gauge boson On the right hand side we already see a factor e −πh µ/2 . It leads to exponential enhancement of mode with h µ < 0, suppress the mode with h µ > 0, and has no effect on the longitudinal mode h = 0. The h µ dependence in the Γ function will be discussed below. Roughly speaking, the mode function receives most of its enhancement when the physical momentum kτ µ where the adiabatic approximation is maximally violated. The produced gauge bosons then follow the comoving dilution in the late-time limit. In Fig. 3 we show this behavior by plotting the mode function A (−1) in (13) with fixed k = 1 against the conformal time kτ . We show the mode functions for different choices of the mass ν and the chemical potential µ. One can see that the early-time behavior as oscillations in kτ with amplitudes and frequency independent of either µ or ν. At late times, the mode functions develop oscillations in log |kτ |, with the frequency determined by ν, and the amplitude determined by both ν and µ. (Note the change of kτ coordinate from linear scale to logarithmical scale at kτ = −10.) The enhancement of the amplitude at later times is evident for large µ.

Approximating the propagators
Later in this paper we will calculate the oscillatory signals from these gauge boson modes. However, it will be helpful first to estimate the signal size before a detailed calculation. For this purpose we use the fact that the signal is contributed mostly by the late-time part of the mode function, and thus will take a look at the late-time expansion of the gauge boson's propagator.
The propagator is constructed following the standard Schwinger-Keldysh (SK) formalism [50], At late times |kτ | 1, we can decompose the propagator into the "local" part and "nonlocal" part, where the latter contains non-integer dependence of the momentum.
Using the asymptotic expression of Γ-function log Γ(z) ∼ (z − 1 2 ) log z − z + 1 2 log 2π, or more explicitly, we can find an expression for non-local propagator for large µ and ν. Note that the asymptotic behavior holds well even with mildly large |b| 1, so we can apply it to, e.g., Γ( 1 2 + i µ − i ν) where the two large numbers µ and ν are cancelling each other. Now we assume µ > 0 without loss of generality, then the h = −1 component is dominant, with the following propagators.
From these results we can read a rule for simple estimate. To summarize this rule more compactly, we introduce the following two quantities, Then the estimate goes as Here we have dropped a prefactor ν −1 in estimating the nonlocal part of D > , since this factor is usually compensated by a positive power ν coming from the in-in integral. In the large mass limit, ν ∼ m A /H, D > local ∝ m −1 A e 2πΩ(γ) . As discussed in App. A, the time integral will change the mass dependence to m −2 A and thus reproduce the EFT limit. At the same time, the exponential factor e 2πΩ(γ) is still present in this limit as the particle production only depends on the relative size of µ and m A . Hence, in the large m A limit, we can approximate a hard propagator as m −2 A e 2πΩ(γ) .

Estimate of Signal Size
With the rule of approximating gauge boson propagator as derived above, we can now estimate the size of gauge boson signals in NG. Different from the fermion case where there is only one relevant 1-loop diagram, here we have at least two relevant categories of diagrams to consider, shown in Fig. 4. Figure 4: One-loop diagrams contributing to gauge boson signals. The blue color marks the soft lines in the squeezed limit.
In the first category shown in Fig. 4(a), the gauge boson loop is directly attached to the inflaton external lines, with the coupling from φF F . In the second category shown in Fig. 4(b) and (c), we mix the Higgs σ (dashed lines) and the inflaton (solid lines) via the dim-6 operator. The mixing is of orderφ 0 σ 0 /Λ 2 Σ . Bothφ and σ 0 could be large so this mixing could be large, too. There are suppressions from Higgs internal lines ∼ 1/m 2 σ but this can well be compensated by the two-point mixing as u → 1. Of course there can be mixed case where the gauge boson loop is attached to both σ and δφ. For now we will only focus on the diagrams in Fig. 4.
We will now estimate the NG of each diagram. We will consider both the non-signal part of the NG, which we refer to as the "background," and the oscillatory signals. The background estimation is useful when translating the current NG constraints to that of model parameters. The way of estimating NG can be roughly summarized as "loop factors × propagators × vertices" multiplied by a numerical factor 1/(2πP 1/2 ζ ) where P ζ 2 × 10 −9 is the amplitude of the scalar power spectrum. We refer readers to [19] for a more detailed description about the signal estimate. We note that there are known subtleties about this estimate. In particular, we can use this rule to figure out easily the exponential dependence on the model parameters (the chemical potential and the mass), but generally we cannot get the correct power dependence. In fact, we don't even know how to calculate this power dependence analytically in general situation. We discuss these subtleties in more detail in App. A. Now we estimate the signal size for each diagram. To simplify our expressions, we will take the unit H = 1 in the rest of this subsection. We will also focus on large mass region ν 1 so that ν m A . Our final results depend only on three free parameters ≡ µ/m A , m A , and u introduced in (3). The diagram in Fig 4(a) can be estimated as We have included a loop factor 1/16π 2 , the prefactor 1/(2πP 1/2 ζ ), three vertices, each of which gives 1/Λ F , and finally the internal gauge boson propagators taken from (23). We estimate two types of NG. f (bg) NL is obtained by approximating each of the internal gauge boson propagators with the local form m −2 A e 2πΩ(γ) . This is the kind of NG which will be constrained by the current observation. f (osc) NL is the size of the oscillatory signal. To estimate this, we approximate the hard line (black in Fig. 4) with the local form propagator, while the soft lines (blue in Fig. 4) should take the non-local propagator in (23).
The sizes of the rest of the diagrams in Fig. 4 can be estimated in a similar way. For Fig. 4(b): Here we have taken the σ-propagators to be 1/m 2 σ since we always focus on m σ 1. Finally, Fig. 4 To get an idea of overall signal strength and also the relative importance of difference diagrams, we can look at a special case where 1 and m A 1. This is the parameter range we are mostly interested in. One might want to consider the case with > 1 where the large exponential factor can lead to great enhancement. However, a large exponential factor could be severely constrained by several physical considerations as we will elaborate in the next section. Therefore we will take those exponential factors as O(1) for now. Then, we see that each diagram is simply a factor of P ζ multiplied by some powers of m A , with an expected range 1 m A φ 1/2 0 60, as well as by a factor of u −3 for Diagrams (b) and (c). The factor P ζ 10 −9 would make the signal tiny, unless there is a large positive power of m A or if u 1. We see that Diagram (a) is independent of u, and the signal f (osc) NL has only one positive power in m A . Therefore Diagram (a) will be tiny in any case. On the other hand, (b) and (c) get more powers of m A . So they will be the dominant contribution when u ∼ 1 or u < 1, although they will also be suppressed when u 1. We note that this result is independent of the relative size between Λ F and Λ Σ as long as all parameters are within the range of validity of our estimate. (In particular, we must have m A > H for our estimate to be valid. Therefore we require σ 0 /H > 1/g.) To summarize, we find that the largest signal is from Diagrams (b) and (c). Even without a truly exponential enhancement, we are able to get large signals with the help of the factor m 4 A , although we should note again that it is generally difficult to estimate to power dependence correctly. Indeed, a more careful calculation in the following section shows that this power dependence is actually m 11/2 A rather than m 4 A . But even this result may not capture the full power dependence on the mass m A . We comment on this issue in the App. A and leave a possible improved calculation for future study.
And it can be observed at this point that the large signal should survive all the constraints mentioned above. This is because all those constraints are actually constraining a large exponential factor, namely that e 2πγ cannot be much greater than 1. But here we see that we do not need a large exponential factor.

Explicit Calculation
The estimate above shows that Diagrams (b) and (c) can possibly give rise to visibly large signals while Diagram (a) is always tiny. Therefore we will calculate (b) and (c) more explicitly in this subsection. We will present some detailes of calculation of Diagram (b). Diagram (c) can be calculated quite similarly, to which we can simply adapt the result of Diagram (b) with appropriate changes. We follow the method in [12] but with improvements. More details about the formalism and the convention we used here are reviewed in [50]. Uninterested readers can directly go to the final results in (47) and (48).
According to the diagrammatic rule, Diagram (b) can be written as Here a i = ±1 (i = 1, 2, 3) are in-in contour indices, and the momenta are as labeled in Fig 4(a). D a i a j µ ν 's are gauge boson propagators, and we have written the external lines compactly in terms of mixed propagators [50], in which G a is the bulk-to-boundary propagator of the inflaton fluctuation and D a± is the propagator for scalar σ. The integrals in (30) are difficult to be carried out directly, and we adopt several approximation to make progress.
Approximation of the mixed propagator. First, since the oscillation signals associated with the mixed propagator are small, we will ignore them. Hence, we can expand the mixed propagator in the large mass limit. This is equivalent to replacing the σ-lines by effective vertices 1/m 2 σ . More explicitly, we derive the effective vertex as follows.
Soft gauge boson propagator. Next, to carry out the loop integral, we will make a late-time (|kτ | 1) expansion to the soft gauge boson propagators (blue lines in Fig. 4). For this purpose, we decompose the gauge boson propagator into components with fixed helicity, and use the nonlocal part of the propagator D > is reviewed in [50]). This is the crudest approximation among all we make in this calculation. There are two unsatisfactory points about it. First, the late-time expansion holds well only when |kτ | 1. On the other hand the time integral in (30) receive contributions for all |kτ | | µ| which is outside the range of validity of the expansion. Second, the local part of the propagator can actually contribute to oscillatory signals, a point that was overlooked in previous studies. We provide a detailed discussion of these issues in App. A.
Hard gauge boson propagator. For the hard loop line (the black line in Fig. 4), [12] evaluates it at the saddle point of the integral while [16] used an "improved" late-time expansion. None of these worked perfectly. So in this work we will work in the large-mass expansion, assuming that this hard line can be approximated by an EFT vertex 1/m 2 A . The advantage of this approximation is that it allows us to carry out the loop momentum integral completely, as opposed to previous works where the loop integral is evaluated at some particular momentum configurations. Therefore, we will make the substitution A µ (x)A ν (y) → g µν e 2πΩ(γ) m −2 A δ (4) (x − y). We have inserted an additional factor e 2πΩ(γ) following (23). This takes account of the fact that the gauge boson mode function becomes tachyonic when γ > 0 and the exponential enhancement applies to both the positive-and the negative-frequency parts of the mode function. Clearly we should not trust this substitution when γ 1. From this we get the effective vertex, Note that we have separated all scale factors so the metric appeared here is the Minkowski metric η µν rather than the FRW metric.
We will also present a less rigorous derivation within the current formalism in App. A.
Simplified Result with Approximations. Above we have listed all approximations we made in order to carry out the integral analytically. As a result, the integral (30) can be recast into the following form where the loop integral I(k 3 , τ 1 , τ 3 ) is defined as Note that we have dropped the SK indices for the soft gauge boson lines. Now this has the same form as Fig. 4(c) after contracting the hard internal line. Working out the vertex coefficients shows δφ 3 Loop integral. Now we perform the loop integral. Without loss of generality we assume µ > 0 in which case the negative helicity component dominates the result. The loop integral can then be written as where p ≡ q − k 3 . Using rotation symmetry we can put k 3 = (0, 0, k 3 ), q = (0, q sin θ, q cos θ), p ≡ (0, p sin χ, p cos χ) = (0, q sin θ, q cos θ − k 3 ). (39) From this we find So the polarization product in the loop integrand is [1 + cos(θ − χ)] 2 /4. Then the integral can be carried out analytically. We keep the non-local part only, and the result is, Time integral. Now it is straightforward to finish the time integral in (35), and the result is It is conventional to represent the result in terms of the dimensionless shape function S(k 1 , k 2 , k 3 ), related to the correlation function δφ 3 through The prime on the correlator means the δ-function of momentum conservation stripped, which was implicitly assumed in all previous expressions. Therefore we have, in the squeezed limit k 1 k 2 k 3 , Similarly, It is useful to show the limit of this shape function when µ ∼ m A H. Using = µ/m A 1 and ν m A /H, we find, This agrees with the previous estimates (27) and (29) in overall parametric dependences, except for the power dependence in m A . The calculation here yields m 11/2 A while the previous estimate gives m 4 A . From the above calculation it is clear that the additional powers m

3/2
A are generated during performing the loop integral and the time integral, which can in no way be easily estimated.

Constraints
Unlike the chemical potential for fermion production studied previously, the gauge boson signal obtained here in (47) and (48) can be exponentially large if naively extrapolated to large 1. Clearly we should not trust the result with a large exponential enhancement as the perturbation theory breaks down there. In addition, with the presence of exponential enhancement, we will also meet several physical constraints. In this section we will outline the constraints in the parameter space.
First of all, there is a constraint from the validity of EFT expansion used in the Lagrangian in Eq. (2). To make sure the derivative expansion in ∂φ is valid when evaluated with the inflation background, we require Λ Σ , Λ F >φ In addition to the above EFT constraint, there are three constraints we will consider. First, the energy density of the produced gauge bosons should be subdominant during the inflation. Second, the equilateral non-Gaussianity should be within the current limit. Third, the perturbative expansion should be justified for our computation of signals to be valid. Basically all these constraints require that the chemical potential not to be greater than the mass of the gauge boson. But there can be slight difference in numerical factors in each case which we shall go over below.
Energy density. One physical constraint is that the gauge boson must not dominate the energy density, ρ A 3M 2 Pl H 2 . That is, we are assuming a conventional inflation scenario where the exponential expansion is driven by the inflaton's potential energy. This is not necessary. In fact one can consider a scenario "warm inflation" from such gauge boson production [27,28,51]. We will leave this for a future study.
The energy density ρ A of produced gauge boson can be found from the mode function in the late-time limit (15). This late-time behavior can be compared with the normalized basis From this we can read the Bogoliubov coefficient β h as Then, the number density n k of produced gauge boson in the phase space is given by Apparently there is no k dependence in this expression, and this is because we are considering the late-time limit where the particles are fully non-relativistic. The real gauge bosons are produced when the physical momentum |kτ | = µ. After the production they are redshifted to all lower wave numbers. So we will restrict the above momentum integral within a sphere of radius µ. Then we see that the gauge boson energy density is We note again that the late-time expansion (15) we used here is valid only when |kτ | 1 while we have performed the k integral of the resulting particle density n k up to k = µ. We expect some corrections to this result when µ H which is probably insignificant. Then we impose the constraint that ρ A is much smaller than the energy density during inflation, One may consider a stronger condition that the gauge boson energy density ρ A should not affect the slow-roll potential of the inflaton. But this is unnecessary since we can always adjust the inflation potential a little bit so that the resulting power spectrum agrees with data even including the effects of gauge bosons.
For not too different from 1, namely m A µ =φ 0 /Λ, we can see that the mass of the gauge boson can never be greater thanφ 1/2 0 by the EFT constraint (49). So m 4 Pl H 2 , where is the first slow-roll parameter. So the density ρ A is never larger than 3M 2 Pl H 2 without the exponential factor. Therefore it is only important to consider the parameters with µ > ν or equivalently > 1. In this case the density is roughly We also note that this bound becomes less constraining for lower scale inflation since we have shown that ρ A 4 e 2π( −1)m A /H ρ by EFT bound. For lower scale inflation both ρ and ∝ ρ get smaller, and thus one could tolerate a larger exponential factor.
Equilateral Non-Gaussianity. As mentioned in Sec. 3.3, the same processes that contribute to the gauge boson signal also generate a "background," which in the EFT limit contribute to the equilateral NG. Equilateral NG has been constrained by Planck measurement to be f (eq) NL O(10). Therefore, using our estimate (26), we have the following constraint, Perturbativity. The gauge boson signals were calculated using perturbation theory. At the diagram level, the gauge boson production manifests itself through the exponential enhancement of gauge boson propagators. A large exponential factor could invalidate the perturbative expansion in terms of diagrams. For example, we can consider the interaction vertex g 2 σ 2 A 2 . Inspection of loop expansion with this vertex shows that the effective expansion parameter is Pertubativity requires (in the limit µ, m A H) Of course this bound depends on the interactions. One could also consider other interactions such as φF F which is less important than the above one. We also note that the breakdown of perturbation expansion is not a physical problem. There could interesting effects in the strongly coupled region of parameter space which we leave for a future work.

Discussions
The main results of this paper are presented in Fig. 1 and 2. We see that in the parameter region of µ ∼ m A H, there is a window of opportunity in which the oscillatory signal of gauge boson production is observable by current or near future probes. f (osc) NL can be as large as O (10). In addition, the gauge boson production will also produce NG in the equilateral limit. Parameter space with even larger oscillatory signal is already constrained by the current limit on the equilateral NG. It is also noteworthy that the signal considered in this paper occupies the oscillation frequency region 4 ν 20, with f (osc) NL 0.1, shown in Fig. 2. This is a signal region distinct from models in the category of quasi-single-field-inflation, and those in which the inflaton has a coupling to the chiral current of massive fermions.
In this paper, we have also developed a set of rules to estimate the size of the signal, with the exponential factor properly taken into account. This helped us to focus on a set of diagrams with dominant contribution. Through a more careful study of the late time behavior of the propagator and their time integral in App. A, we also notice that the part of the propagator which is analytic in momentum can also contribute to the oscillatory signal, which has missed by previous studies. We also found that while the exponential factor can be reliably estimated, the power dependence on the mass parameters (gauge boson mass m A in our case) can not be captured by a simple counting. While the approximate calculation in Section 3.4 improves upon the simple estimate, it is unlikely to contain the fully accurate power dependence on m A . This underscores the important to perform a full calculation without approximation in order to set completely accurate constrains and make projections. This is a promising direction to pursue in the future.
Acknowledgment. We thank Xingang Chen, Junwu Huang, Matt Reece, and Yi Wang for useful discussions. LTW is supported by the DOE grant DE-SC0013642.

A More on Estimating the Signal Strength
In [19] we showed how to estimate the signal strength by assigning simple factors to each vertex and propagator in the SK diagram. This simple shortcut can get correctly the exponential dependence from propagators as well as the power dependence from vertices. However, it is generally difficult to get the power dependence from the propagator correctly. In this appendix we elaborate on this issue. The main conclusions are: 1. It is not sufficient to include power dependence on the mass parameters from the propagators and the vertices, since the SK time integrals can introduce additional power dependence (although these integrals do not introduce additional exponential dependence.). It seems there is no simple way to count the powers from the time integrals.
2. For heavy propagator with large mass m H that can be approximated as EFT local operators, one can show that the power dependence on the mass is always 1/m 2 , by either EFT argument or explicit calculation. There could also be exponential dependence for EFT part due to chemical potential which can also be reliably estimated.
3. For non-local "on-shell" propagators, there is no simple way to estimate the power dependence correctly. It is likely that even more careful calculation with late-time expansion of propagators cannot get the powers right, either. However, we expect the late-time expansion can at least capture a fraction of oscillation signals, so it is still a useful way to evaluate the signal strength.
To illustrate these point we will first consider a simpler case in flat space, and then go to the inflation background.
Mass dependence in flat-space correlators. First we will show that the time integral can generate additional power dependence. For this purpose it is helpful first to look at the flat-space example. In flat space there is no oscillation features in the correlation function, and thus we will only consider the "local" part.
In flat space the propagator of a scalar field σ of mass m σ can be written as where E k = m 2 σ + k 2 . In the large m limit the propagator goes like 1/m, the same as the propagator in inflation. So naively we would estimate the contribution of such a propagator as 1/m for large m. However, we know that the correct EFT limit is 1/m 2 . So where is the additional power dependence from? 3 The answer is that the time integral will contribute another power. To see this, we calculate the correlator of 4 light scalar fields φ at t = 0 connected by an s-channel σ, with m φ m σ and vertex λφ 2 σ. Using the diagrammatic rule, the correlator is Here E i = k 2 i + m 2 φ (i = 1, · · · , 4) and E s = k 2 s + m 2 σ . Here we can already see that the time integral will introduce additional power dependence on m σ at large m σ such as dt 1 e iEst 1 ∼ E −1 s ∼ m −1 σ . However, superficially we would expect two powers of m −1 σ being introduced since we have two time integrals. This combined with m −1 σ would give 1/m 3 σ behavior, not in agreement with EFT counting 1/m 2 σ . But if we look at the integral more closely, by writing I = λ 2 (I T + I N ) in terms of the time-ordered part I T and non-time-ordered part, I N , we will get, where E 12 ≡ E 1 + E 2 and E 34 = E 3 + E 4 . So it is the time-ordered part I T ∼ 1/m 2 σ that gives the correct EFT behavior for large m σ , while the non-time-ordered part I N ∼ 1/m 3 σ gives subleading contribution. Our previous naive guess implicitly assumed that the two time integrals can be factorized and thus applies only to I N .
Mass dependence in inflationary correlators. The above example shows that the timeordered integral is important to get the correct power dependence. Similar calculation in the inflation background can also be done, giving similar result. As a demonstration, we use the example of quasi-single-field inflation with one massive field σ. The relevant operator is (∂φ) 2 σ. Expanding around background value, φ = φ 0 + δφ and σ = σ 0 + δσ, it gives rise to a 3-point vertex (δφ ) 2 δσ and a two point mixing δφ δσ. Putting these together, we can form a 3-point diagram with one massive propagator σ in the middle. From naive estimate in the large mass limit ν = (m 2 σ /H 2 − 9/4) 1/2 1, the leading EFT piece scales like ν −2 , while the signal scales as e −π ν .
Analogous to the flat space example, we will separate the integral into a time-ordered piece 2. The result should be compared with the full result in [52], in which the full propagator was used without making late-time expansion. Then in the large mass limit ν 1, the signal scale as ν 3/2 e −π ν instead of the late-time result νe −π ν . The mismatch of the powers signifies a failure of the late-time expansion.
3. The second term in (68) shows that the local part of the propagator (which is analytic in k) can also contribute to the oscillatory signal after time integral. A similar result was also observed in [52] for the more suppressed signal in the non-time-ordered integral.