The Scalar Chemical Potential in Cosmological Collider Physics

Non-analyticity in co-moving momenta within the non-Gaussian bispectrum is a distinctive sign of on-shell particle production during inflation, presenting a unique opportunity for the"direct detection"of particles with masses as large as the inflationary Hubble scale ($H$). However, the strength of such non-analyticity ordinarily drops exponentially by a Boltzmann-like factor as masses exceed $H$. In this paper, we study an exception provided by a dimension-5 derivative coupling of the inflaton to heavy-particle currents, applying it specifically to the case of two real scalars. The operator has a"chemical potential"form, which harnesses the large kinetic energy scale of the inflaton, $\dot{\phi}_{0}^{1/2} \approx 60H$, to act as an efficient source of scalar particle production. Derivative couplings of inflaton ensure radiative stability of the slow-roll potential, which in turn maintains (approximate) scale-invariance of the inflationary correlations. We show that a signal not suffering Boltzmann suppression can be obtained in the bispectrum with strength $f_{\mathrm{NL}} \sim \mathcal{O}(0.01-10)$ for an extended range of scalar masses, $M \lesssim \dot{\phi}_{0}^{1/2}$, potentially as high as $10^{15}$ GeV, within the sensitivity of upcoming LSS and more futuristic 21-cm experiments. The mechanism does not invoke any particular fine-tuning of parameters or breakdown of perturbation-theoretic control. The leading contribution appears at tree-level, which makes the calculation analytically tractable and removes the loop-suppression as compared to earlier chemical potential studies of non-zero spins. The steady particle production allows us to infer the effective mass of the heavy particles from the variation in bispectrum oscillations as a function of co-moving momenta. Our analysis sets the stage for generalization to heavy bosons with non-zero spin.


Introduction
The paradigm of cosmic inflation (for a review see [1]) gives a robust mechanism for the high degree of homogeneity and isotropy of the universe on very large scales, as the result of exponential spacetime expansion driven by classical inflaton scalar dynamics coupled to General Relativity. Furthermore, quantum fluctuations in this paradigm lead to a structure of tiny inhomogeneities, beautifully consistent with what is observed, for example, in the Cosmic Microwave Background (CMB). In particular, data suggests that these primordial fluctuations are scale-invariant, adiabatic and Gaussian to a good precision [2,3], supporting their origins from a single weakly coupled quantum field, minimally identified with the inflaton itself. In general, interactions can give rise to small non-trivial n-point correlations of the inflaton field, which translates into non-Gaussianity (NG) of the primordial curvature fluctuations (for reviews see [4,5]). But with the expected improvement in experimental precision in the near future, especially with Large-Scale Structure (LSS) data [6,7], we can hope to detect small NG, revealing subdominant interactions of the inflaton with itself or other fields. In this way, the properties of these other fields during inflation may be encoded in primordial NG observables. 1 In more detail, the time-dependent inflationary space-time provides energy of order the expansion rate, i.e., the Hubble scale (H), that leads to the on-shell production of any particle with mass M ∼ H. The particle production here is associated with the fact that a geodesic observer in a de Sitter (dS) space-time sees a thermal distribution of particles at the Hawking temperature T Hawking = H/2π (see e.g. [8]). If these particles interact with the inflaton, once produced they can later decay into inflaton fluctuations, resulting in primordial NG [9]. The cosmological production and later decay into inflatons of heavy particles leaves a unique non-local feature in the primordial 3-point function of curvature fluctuations (R), also known as the bispectrum. When expressed in co-moving momentum space, this spacetime non-locality appears as momentum non-analyticity in the so-called squeezed limit, where one of the momenta is soft, k 1 ∼ k 2 k 3 . Typically [9][10][11], (1.1) Here s is the spin of the heavy particle and θ is the angle between the hard and the soft momenta. We see that the non-analytic dependence on the momentum ratio k 1 /k 3 , which is either oscillatory (M > 3H/2) or non-standard power-law (M 3H/2), and the angular dependence can be used to do spectroscopy of the mass and the spin of the heavy particle. This opens up a new window into particle physics at energy scales as high as the inflationary Hubble scale ( 5 × 10 13 GeV [2]), potentially orders of magnitude beyond the reach of any terrestrial collider. This program of research and its applications have been dubbed "cosmological collider physics" . While this is a very exciting program, a significant hurdle to fully exploit the reach of the "cosmological collider" is that the non-analytic signal gets exponentially smaller (∼ e −πM/H ) for M H as seen in eq. (1.1). Intuitively, this is a "Boltzmann suppression" factor due to the thermal-like nature of particle production at the Hawking temperature T Hawking = H/2π, where the bispectrum gives the amplitude for particle production, ∝ e −M/2T Hawking , corresponding to the standard Boltzmann probability e −M/T Hawking . On the other hand, for small masses M H, we see that the non-analytic factor in k 1 /k 3 becomes approximately analytic and difficult to disentangle from other NG sources, such as inflaton self-interactions. Therefore apparently, these two considerations very strongly restrict observable masses to a window of ∼ H.
Of course, in addition to the production amplitude for new heavy particles, the strength of NG depends on the couplings to the inflaton, which we now consider. Radiative stability of the inflaton potential strongly suggests that it should have derivative couplings predominantly, schematically of the form where O heavy is made from the heavy fields and Λ is the scale at which the non-renormalizable EFT description breaks down. While the exact cutoff depends on the assumptions about the nature of the UV-physics, plausibly it should be at least as big as the scale of inflaton kinetic energy in slow-roll inflation [44], where φ 0 (t) denotes the classical inflaton field trajectory, and where we have used the fact that the CMB temperature power spectrum amplitude implies H 2 /φ 2 0 ≈ 10 −7 [2]. It therefore seems that the typical (dimensionless) couplings cannot offset Boltzmann suppression, since they are at best O(1) when some of the inflaton fields are evaluated at the classical expectationsφ 0 , or are significantly suppressed.
The case n = 1 in eq. (1.2) is exceptional however, with the inflaton coupling to a current made of heavy fields, ∂ µ φJ µ heavy /Λ. This has been studied in the context of heavy fermionic matter [26,33,34,37,45,46] and heavy spin-1 bosons [37,41,47,48]. Here we will study the simple but interesting case of heavy scalar matter, and show that there are dramatic new features enhancing the mass reach of cosmological collider physics in a controlled and robust manner. The case of a single real heavy scalar, σ, is however trivial since the unique coupling is given by Integrating by parts and replacing the resulting φ by −V (φ) using the inflaton equation of motion, this only makes a mass contribution to σ modulo small slow-roll corrections. For the case of two real scalars, σ 1 , σ 2 , on the other hand, eq. (1.4) generalizes to For generic values of Λ 1 and Λ 2 , these operators can not be eliminated by integration by parts. Inserting the time-dependent inflaton VEV generates a quadratic mixing between the two scalars, where λ 1 =φ 0 /Λ 1 and λ 2 =φ 0 /Λ 2 are constant in time up to slow-roll corrections. This form is reminiscent of the chemical potential for a complex scalar χ = σ 1 + iσ 2 . Indeed, for the special case λ 2 = λ 1 = λ and equal masses M 1 = M 2 , eq. (1.6) is the chemical potential ∼ λJ 0 , where J is the Noether current associated to χ phase rotations, J 0 =σ 1 σ 2 −σ 2 σ 1 . We know that the presence of a chemical potential in standard thermal equilibrium can lead to occupation of excited states with energy greater than the temperature, thereby overcoming Boltzmann suppression. Applying this intuition to the thermal-like nature of particle production in dS spacetime, we expect particles with masses M λ to be produced efficiently. 2 Let us estimate the value of chemical potential in our case. From eq. (1.3), we see that λ φ 0 /Λ φ 1/2 0 ≈ 60H, which is T Hawking . That is, we can expect that particles much heavier than H can be produced without the dS analog of Boltzmann suppression, greatly widening the window of observable masses in the cosmological collider physics program.
In this paper, we demonstrate that the above expectations of non-Boltzmann-suppressed NG mediated by heavy scalars H M λ can indeed be realized. The simplest way to see how this works is to first note that the "chemical potential" coupling can be thought of as χ being charged under a "gauge" field which is pure gauge in form, A µ ≡ ∂ µ φ/Λ. It can therefore be removed by doing a gauge transformation on χ. However, this will induce gauge transform phase factors in any other gauge-non-invariant masses or couplings in the theory: where O is an interaction/mass in χ violating the fake "gauge invariance" by q units (that is, there are q more powers of χ than χ † ). Given the slow-roll approximation, φ 0 ≈φ 0 t, we see that we get effective high-frequency couplings of the heavy fields.
For example, if we depart from our assumption of equal heavy masses, or from λ 2 = λ 1 , there are effectively χχ-type symmetry-breaking mass corrections accompanying the chemical potential, which will be multiplied by e 2iλt after the field redefinition. In general, in this way we get time-dependent χ couplings and mass terms with frequencies of order λ. These couplings can then naturally result in efficient particle production for masses up to λ H. 3 We will study the case of O(χ, χ † ) with q = 1 and show that it can lead to striking NG signatures of the heavy scalars, the largest and the most theoretically tractable occurring via tree-level exchanges with the inflaton. These can be readily observable in upcoming LSS experiments [6,[49][50][51] and future 21-cm cosmology [52,53]. We will also briefly study the case q = 2, which is mostly closely analogous to the fermionic case. In general, both q = 1 and q = 2 can co-exist, but for simplicity, we pursue them separately.
There are some qualitative differences between the heavy scalar case here and the earlier studies of heavy fields with non-zero spin s. In these previous studies [26,33,34,37,37,41,45,47,48], the chemical potential modifies the dispersion relation with a term linear in momentum using a coupling ∝ k · s. The modified frequency ω(t) violates the adiabaticity condition, i.e.ω/ω 2 1, at a certain time during the evolution of mode functions, leading to particle production. Of course a similar effect is impossible for spin-0 fields. While these spin dependent effects are absent for scalar fields, a chemical potential can still manifest as an explicit time-dependence in mass corrections and couplings as seen in eq. (1.8), which is the source of particle production. Ref. [37] surveyed different chemical potential possibilities, but concluded that there was no effect for complex scalars as given by eq. (1.5). In this paper, we will show that there is indeed a strong effect within theoretical control.
Our main focus will be on q = 1 case with O = χ, which has the advantage that the contribution of the inflaton-scalar coupling to the bispectrum occurs at tree-level. This removes the loop-suppression seen in gauge boson and fermion cases, gives a robustly large mass reach M φ 1/2 0 ≈ 60H, while allowing an analytically tractable calculation. Note this is significantly higher than the chemical potential reach for fermions φ 1/4 0 H 1/2 ≈ 10H. While the spin-1 chemical potential reach is as high as the scalar case M φ 1/2 0 ≈ 60H, getting a viable signal requires a period of exponential growth in particle production [47], in turn needing a very close coincidence of the mass and chemical potential scales in order to stop in time. By contrast in the scalar case, no exponential growth or coincidences are necessary.
The field redefinition leading to the form of Lagrangian, eq. (1.7), demonstrates that the scalar chemical potential is equivalent to a special form of particle production due to non-derivative couplings to the inflaton, similar in spirit to [54]. However there are some significant differences. Generally, non-derivative couplings will induce radiative corrections to the inflationary effective potential upon integrating out the heavy matter, of concern given the fragility of slow-roll inflation. Furthermore, non-derivative couplings of (nearly) massless scalars ordinarily yield divergent late-time correlators. However, in our case the non-derivative couplings to the inflaton do not break the shift symmetry under which φ → φ+c when it is accompanied by the phase rotation χ → e −ic/Λ χ. This means that at the level of the inflationary effective potential, radiative corrections from integrating out matter must respect φ shift symmetry, and hence vanish. This matches the obvious conclusion from the original form of interaction, eq. (1.5), where the derivative couplings of φ cannot generate an effective φ potential. Similarly, the manifest derivative couplings in eq. (1.5) ensure that all late-time divergences cancel when the correlators are computed using eq. (1.7) (which however is more useful in computing the effects of the chemical potential). Finally, it may appear worrying that the non-derivative couplings in eq. (1.7) can be Taylor expanded in powers of φ/Λ even though the φ field transit is much greater than Λ. In general, nonderivative couplings containing (φ/Λ) n would not have a controlled expansion. But here the symmetry requirement that the inflaton factor non-linearly carries "charge" q uniquely gives the re-summed form e iqφ/Λ , thus matching the manifest EFT control in eq. (1.2).
As seen in eq. (1.8), the particle production resulting from the chemical potential is continuous and approximately constant during inflation, given by a steady frequency of the effective couplings/masses, as opposed to a more time-dependent or punctuated particle production. Relatedly, the shift symmetry gives us an approximate time translation invariance, which combined with dS isometries implies that approximate scale invariance of our correlators is maintained. This structure of correlators is useful in doing the spectroscopy of the heavy particles from the bispectrum in a simple way.
The rest of the paper is organized as follows: we describe the relevant observables and fix the notation in section 2. Then in section 3, we present our model that illustrates the mechanism of chemical potential in the context of heavy complex scalars. Section 4 contains the main results of this paper where we calculate the bispectrum in the squeezed limit and analyze the strength of NG in various regions of the parameter space. There we also give a simple and intuitive estimate of the bispectrum using the method of stationary phase while a full calculation is carried out in appendix B using the results from appendix A. After discussing in section 5 a novel procedure to infer the effective mass of the heavy particle during inflation from the momentum dependence of the bispectrum, we conclude in section 6.

Preliminaries
Let us start with the definitions and the notation that we will use throughout the paper. We will work in (−, +, +, +) signature of the metric. Under the slow-roll approximation, the fractional change in the Hubble constant |Ḣ/H 2 | is 1% [2]. Since H changes slower than other quantities during inflation, we take it to be a constant and approximate the metric during inflation by the dS metric, where η is the conformal time. It is related to the proper time t through the relation dη = dt/a(t), where a(t) is the scale factor. This gives ηH = −e −Ht during inflation. We set H = 1 in the rest of the text for the ease of calculation, unless it is explicitly written.
Factors of H can be restored by dimensional analysis. The inflaton field can be separated into a classical homogeneous background φ 0 that sources inflation, and the (quantum) fluctuations ξ that seed density perturbations, φ = φ 0 (t) + ξ(t, x). In the spatially flat gauge [55], the fluctuations ξ(t, x) are related to the gauge-invariant comoving curvature perturbation R (for a review and original references see e.g. [56]) as The observations require the primordial power spectrum P R (k) to be (almost) scale invariant, where n s ≈ 0.96, k * = 0.05 Mpc −1 is the pivot scale, andφ 1/2 0 ≈ 60H [2]. Here and below, the denotes the convention that momentum conserving factors have been taken out, R k 1 · · · R kn ≡ (2π) 3 δ 3 ( k 1 + · · · + k n ) R k 1 · · · R kn . (2.5) The primordial non-gaussianity (NG) corresponds to having connected non-zero n−point correlation functions of R. More specifically, the bispectrum is the 3-point function in momentum space denoted by B(k 1 , k 2 , k 3 ), It is a convention to denote the bispectrum by F (k 1 , k 2 , k 3 ), which can be thought of as the bispectrum appropriately normalized by the power spectrum . (2.7) In our analysis, we will be mostly interested in the squeezed limit, i.e., when k 1 ≈ k 2 k 3 . Using eq. (2.3) and eq. (2.4), the expression for F can be simplified in this limit as It is customary to typify the strength of NG by a single number, f NL , at the "equilateral" configuration where all the momenta have the same magnitude, This definition is consistent with the convention for local NG (f local NL ) at the equilateral point. Using eq. (2.3) and eq. (2.4), this can be written in terms of inflaton correlation functions, where | k 1,2,3 | = | k|.
The constraints on f NL vary depending on the shape of the bispectrum. The latest Planck analysis sets the bounds as |f NL | ∼ O(5 − 50) depending on the shape of NG [3]. In the near future, LSS observations (see e.g. [6]) are expected to probe δf NL ∼ O(1) -bringing several of our benchmark scenarios under the LSS reach, as we will see later. More futuristic cosmic variance limited 21-cm cosmology experiments can reach even better sensitivities with δf NL ∼ O(10 −2 ) or even smaller [52,57].
We use the in-in formalism to evaluate the bispectrum. The reader may refer to [58] for a more comprehensive review. Using this formalism, the expectation value of a Heisenberg operatorÔ H at any single time t can be perturbatively computed using, where we can expand in powers of the interaction Hamiltonian H I int to any desired order. Here we see the appearance of both time-and anti-time ordered exponentials, respectively denoted by T andT operators. The state |0 is the free Bunch-Davies vacuum but the i prescription projects it onto the interacting vacuum |Ω . The superscript I denotes fields in the interaction picture. The bispectrum is the 3-point correlation function evaluated at the end of inflation, i.e.,Ô(t) → (ξ k 1 ξ k 2 ξ k 3 )| η≈0 in eq. (2.11).

Chemical potential for a complex scalar with U(1) symmetry breaking
Before getting into our actual models, let us consider a simple case of a complex scalar field χ charged under a global U (1), where we introduce an explicit chemical potential by hand. The conserved current J µ = −i(χ∂ µ χ † − χ † ∂ µ χ) can be used to introduce a non-zero chemical potential, λ, by shifting the Hamiltonian density as H → H − λJ 0 . Integrating out the canonical momentum in the partition function, Note that adding a chemical potential in the Hamiltonian is equivalent to introducing coupling of the complex field to an external gauge potential A µ with a non-zero time component A 0 . Then A 0 (≡ λ) is the chemical potential.
In the standard thermodynamic setting, the subsystem under consideration can gain or lose charges to the heat bath, which is key to chemical potential having physical effect. In the inflationary setting on the other hand, the subsystem is the entirety of space and hence, the charge conservation is exact. This can be easily seen by noting that any non-trivial effect of the chemical potential can be eliminated by a suitable field rotation, χ → e −iλt χ. 4 Nevertheless, non-conservation due to the heat bath can be replaced here by small explicit breaking of the U (1) symmetry. Consider the simplest case of adding a linear symmetrybreaking term as follows, After the field rotation as before χ → e −iλt χ, we see that the symmetry-breaking term retains the effect of chemical potential in the form of a time-dependent phase, For two independent real scalar fields (σ 1 , σ 2 ), such symmetry-breaking terms are naturally present when written in terms of χ = σ 1 + iσ 2 .
In our model, we will realize just such a residual time-dependence with frequency λ, capable of producing particles in association with inflatons at energies of order λ. The linear term will source a time-dependent spatially homogeneous vacuum expectation value (VEV) for the complex scalar, which means that the leading NG will be at tree-level. The advantage over earlier discussions of NG with chemical potential is two-fold: we no longer pay the 1/(16π 2 ) loop suppression, and the calculation is a lot more tractable as compared to the loop-level calculations in scenarios involving chemical potential for fermions and gauge bosons.

Realizing a chemical potential via inflaton couplings
Coming to the case of inflation, the chemical potential can be implemented with derivative coupling of inflaton with the complex scalar current J µ . Consider the Lagrangian, Here V (φ) is some suitable inflationary potential whose details we do not need to specify. The second line describes the interactions, out of which the first term is the coupling of the inflaton to the current. In the slow-roll approximation, the inflaton gets a VEV such thaṫ φ 0 (60H) 2 is a constant to a good approximation. This generates a chemical potential, where λ ≡φ 0 /Λ. It is evident from the discussion above that we must also include a small explicit symmetry-breaking term, which is the last term on the second line in eq. (3.4).
The second term, with the coefficient c, should be naturally present as it is consistent with the shift symmetry of inflaton as well as the U (1) symmetry of the complex scalar. The coefficient of this term gets modified as it is also generated from the kinetic term for χ after the field redefinition that removes the dim-5 inflaton-current coupling, Then eq. (3.4) can be written as In the new field variableχ, the chemical potential gets translated into a non-trivial phase in the symmetry-breaking term. Another advantage of this step is that the Lagrangian in eq. (3.7) no longer has interactions with time-derivatives ∂ 0χ . It is then straightforward to get the interaction Hamiltonian starting from the Lagrangian density using H int = −L int . 5 For simplicity, we restrict the coefficient c ≤ 1. Note that c = 1 is a special case that corresponds to the addition of chemical potential without affecting the mass, as seen in eq. (3.1). Here, we keep this coefficient arbitrary.

Analysis of the spatially homogeneous VEV of the complex scalar
In the following, we will restrict our attention to the parameter space where interactions in eq. (3.7) are perturbative such that the inflaton dynamics is not altered at leading order in perturbation theory. Therefore we can treat the homogeneous part of the inflaton as φ 0 whereφ 0 ≈ (60H) 2 . This perturbativity is ensured by taking Λ(α) sufficiently large(small), as discussed in more detail in section 4.2. The term linear inχ in eq. (3.7) breaks the U (1) symmetry explicitly and gives a time-dependent VEVχ 0 which we estimate now. From eq. (3.7), the quadratic Lagrangian forχ is The bare mass M gets modified by the chemical potential to give an effective mass The parameter c controls the contribution of chemical potential to the effective mass. As mentioned before, c = 1 is a special case where the coupling of inflaton to the current only generates a chemical potential and does not affect the mass, i.e., M eff = M . Due to the time-dependent phase in the term linear inχ, the standard procedure of shifting the field by a constant value will not get rid of this tadpole completely. We instead solve the equation of motion (EOM) ofχ 0 following from eq. (3.8). The homogeneous part of the inflaton field during slow-roll can be written as sinceφ 0 (t) constant in the slow-roll approximation. The conformal time η i depends on the initial value of φ 0 . Since η i will only appear as an overall phase factor, its exact value will not affect our results. From eq. (3.8), the EOM for the complex scalar field (in the dS limit of inflationary spacetime) is This is a standard EOM for a massive scalar field in the inflationary background with an extra time-dependent piece on the RHS coming from the linear term. This extra piece will source spatially homogeneous part, i.e., k = 0 mode. We writeχ →χ 0 (η) + δχ( x, η), whereχ 0 is the time-dependent VEV, and δχ( x, η) is the fluctuation. Putting ∇ 2χ 0 = 0 in eq. (3.11), we get the equation for the homogeneous part, The general solution can be written as, where µ = M 2 eff − 9/4. The first term corresponds to oscillations of the massive field given by its mass ∼ µ, and the usual 1/ √ volume dilution of mode functions. Due to the dilution, this term quickly becomes negligible. On the other hand, the second term exhibits a forced-oscillator type behaviour coming from theχ-tadpole in eq. (3.8), which gives time-dependence to the VEVχ 0 . This term does not dilute (nor diverge) in late time and vanishes at the end of inflation asφ 0 → 0. Thus, the homogeneous part of the complex field can be written asχ Notice that the amplitude of the oscillating VEV in eq. (3.14) depends on the detuning between the natural frequency M eff and the external frequency λ, just like the case of a forced oscillator. The damping comes from the Hubble friction ∼ 1 η ∂ ηχ in eq. (3.11).

Quantization and mode functions of the complex scalar
After doing the analysis for the homogeneous solution, let us move on to discuss the fluctuations. Taking a Fourier transform in eq. (3.11), we get the equation for the k-th mode ( k = 0) as Here, we have ignored a quadratic mixing between δχ and ξ. We will later take this into account perturbatively in the insertion approximation. The RHS in eq. (3.11) is purely spatially homogeneous and does not affect k = 0. The equation for δχ( k, η) is the standard EOM for a massive scalar field in the dS limit of inflationary spacetime. In terms of the redefined fieldχ, the effect of the chemical potential is solely on the time-dependent VEṼ χ 0 , while the fluctuations follow the same dynamics as a free massive scalar field. The general solution is Here, H (1 or 2) iµ are the Hankel functions. The quantized field can be taken to have the following form, are the destruction (creation) operators for particle and antiparticle excitations, respectively, i.e. a k |0 = b k |0 = 0, where |0 is the Bunch-Davies (BD) vacuum. Imposing canonical quantization, using the conjugate momentum π( x, η), given that the mode functions satisfy the Wronskian conditionḡ where s denote derivatives with respect to conformal time η.
Treating all the interactions perturbatively, the quadratic part of the Lagrangian for fluctuations corresponds to a free complex scalar in dS space-time. Thus, minimization of Hamiltonian requires choosing appropriate positive frequency mode at early time, i.e., f k , g k ∝ e ikη at η → −∞. These initial condition and the Wronskian condition fix the form of the mode functions g k , f k to bē 19) where N f = √ π 2 e πµ/2 e −iπ/4 and N * g = √ π 2 e −πµ/2 e +iπ/4 . It can be checked that as required by imposing BD initial conditions. Assuming the validity of perturbative analysis (discussed in section 4.2), we can then take the standard M → 0 limit for inflaton. Quantizing inflaton fluctuation as we get the mode functions by setting µ = 3i/2 in eq. (3.19), (3.23)

Boltzmann-like suppression for fields with no chemical potential
Having studied the mode functions, let us take a detour to understand the origin of the Boltzmann suppression in the absence of chemical potential. Consider a heavy scalar field σ with mass M in dS space-time. The free field EOM is the same as in eq. (3.15). The scalar field is then quantized using BD initial condition as, 24) The BD initial condition means that at early time, the modes were deep into the horizon and hence look like Minkowski modes in the appropriate time variable. This is reflected in the early time limit of the mode functions, where the positive frequency solution is chosen as it minimizes the Hamiltonian at early times. In a time-dependent space-time, the same mode may not be the local positive frequency mode at late times. This point becomes apparent in the late time limit where (up to some unimportant phases) where in the second line, we have used Notice that a small negative frequency part, (−η) −iµ ∼ e iµt , has developed, with a relative amplitude e −πµ , which corresponds to 'cosmological particle production'. Then the probability of particle production is ∼ e −2πµ , which is approximately e −M/T Hawking , where the cosmological Hawking Temperature is given by T Hawking = H/(2π). The non-analytic signature in the bispectrum comes from the 'on-shell' production and propagation of heavy particles, which is always accompanied by this exponential suppression unless there is another energy source of particle production. In our model, we create just such a source by harnessing the kinetic energy of the inflaton through its coupling to the complex current. This introduces an additional explicit time-dependence with frequency λ H, and hence produces heavy fields more efficiently.

Interactions of the complex scalar with the inflaton
We focus on the interactions contributing to the bispectrum at tree-level, which is the dominant contribution in our model. These include terms that mix the complex scalar field with inflaton fluctuations, H mix , and 3-point interactions, H 3 . The relevant terms in the Lagrangian can be obtained from eq. (3.7) after expanding the phase ofχ, and using the VEV ofχ from eq. (3.14). The quadratic mixing terms are where The relevant 3-point interaction terms are given by, where Here we have read off the interaction terms directly from the Lagrangian using H int = − L int , which holds at third order in fluctuations (as explained after eq. (3 .7)). The extra time dependence (−η) ±iλ in the vertices is the most important consequence of the chemical potential. It effectively injects (removes) energy of order λ at each vertex, making the production of fields as heavy as M eff ∼ λ possible without the dS Boltzmann suppression.

Approximate calculation of the bispectrum in the squeezed limit
Having worked out the interaction terms, we now evaluate the bispectrum following eq. (2.11), where ξ 0, k is the inflaton fluctuation with comoving momentum k at η → 0. By expanding the exponential operator, it is easy to see that the bispectrum vanishes at zeroth order and the leading contribution comes from tree diagrams. While in general, the bispectrum contains both analytic and the non-analytic pieces, it is the latter that encode distinctive onshell effects of the massive particle. Hence, for the purposes of this paper, we will compute only the non-analytic parts, keeping in mind that for a complete analysis one would also need to take into account the analytic contributions. However, the analytic contributions can be estimated to be comparable to our non-analytic estimates below in the regime of mild squeezing, i.e., analytic and non-analytic contributions to f NL are comparable. The general form can be written as, 2) The '+' indicates that the vertex that comes from e −i H int dη with time ordering, while '−' comes from e +i H int dη with anti-time ordering. Note that the first sign on I ±± corresponds to the cubic interaction vertex, and the second sign corresponds to the inflaton-scalar mixing vertex. The function h(η 1 , η 2 ) captures the explicit time-dependence of the form (−η) ±iλ in the vertex coefficients in eqs. (3.29) and (3.31). The functions d ± ,d ± are the inflaton bulk-boundary propagators, and they can be evaluated using eq. (3.23) as , and k 12 = k 1 + k 2 . G ±± and F ±± are bulk-bulk propagators for the complex scalar field. These can be evaluated using eq. (3.17) as (4.4) The remaining propagators can be related to these using relations, (4.5) We will fix the proportionality factors in the above by using the relevant coupling parameters in a moment. I −+ and I −− are complex conjugates of I +− and I ++ , respectively. Thus, we only evaluate I +− and I ++ explicitly.
We compute the full analytic form of these contributions in appendix B, however next we show the same calculation using the stationary phase method, which is more transparent in terms of demonstrating particle production through chemical potential. It also proves to be a good approximation in the parameter regime of our interest, λ M > H, as we will see in fig. 2.

Dominant contribution: I +− diagrams
Let us start with the evaluation of I +− diagram. For the purpose of illustrating the mechanism, we focus on the sub-diagram consisting of vertices with larger coefficients, i.e., β 2 and ρ 2 . We will see later that the NG estimate thus obtained matches the full analytic calculation very well.
Using the form of vertices in eq. (3.29), eq. (3.31), and eq. (4.2) the sub-diagram has the form (4.6) Looking at eq. (4.4), we see that the first term on the second line corresponds to the contraction δχδχ † , describing particle propagation, while the second corresponds to δχ † δχ , describing anti-particle propagation. Without loss of generality, let is take λ > 0 in the following analysis. We will see that with this choice, the contribution from the term with F +− propagator dominates. It will become clear that if λ < 0, the analysis for F +− and G +− is simply switched, but our conclusions remain unchanged. Let us start with the second term with F +− propagator. There are two integrals, the one with soft inflaton leg (k 3 ) coming from anti-time ordering, while the other with hard momenta (k 1/2 ) coming from time ordering.
where I is the η 1 -integral from the time-ordered part with cubic interaction. Using mode functions for the inflaton and complex scalar from eq. (3.19), eq. (3.23), we get .

(4.8)
For fields with mass greater than the Hubble scale, the rate of expansion of space-time is much smaller as compared to their natural frequency, i.e.,ω/ω 2 1. Thus, the dynamics of the mode functions can be effectively captured by the adiabatic modes. For the ease of calculation, let us make the substitution x = −k 3 η 2 . Then the full mode function can be approximated by the adiabatic mode as, The integrand then separates into a slowly-varying polynomial function and a highly oscillatory phase g(x), where Here, x 0 corresponds to the time at the start of inflation. The dominant contribution comes from the stationary points where g (x * ) = 0, i.e., − λ x * + 1 + 1 + µ 2 x 2 * = 0. In this case, there is only one stationary point at x * = (λ 2 − µ 2 )/(2λ). Using this, we get g (x * ) = 4λ 3 /(λ 4 − µ 4 ) for λ > µ. and the integration can be evaluated as Here we have used the fact that x * 1, true for most of the parameter space we are interested in. The phase factor δ(λ, µ) = g(x * ) − π/2 is independent of k 3 . The additional phase coming from the choice of initial point x 0 is arbitrary, but it cancels with the phase from I (+) k 12 integral. The role of the chemical potential is now apparent. The existence of stationary phase indicates on-shell production of particles. The phase (−η) −iλ coming from the chemical potential can be thought of as providing energy ∼ λ into the vertex, which is converted to produce inflaton and complex scalar fluctuations at time x * . If λ < µ, there is no stationary point and we will get an exponential suppression ∼ e −π(µ−λ) , which as λ → 0, becomes the usual Boltzmann-like suppression e −πµ for heavy fields.
Doing a similar calculation for I (+) k 12 yields, where (4.14) Stationary phase occurs when λ x * − p − 1 + µ 2 x 2 * = 0. Then in the large p limit, for λ > µ, x * ≈ λ−µ p , and h (x * ) ≈ − p 2 λ−µ . This gives, . Since x * depends on p, we get a p-dependent phase from e ih(x * ) , (4. 16) In the above, we have taken a squeezed limit p λ/µ ∼ 1, where this expression can be simplified to separate out the momentum-dependent phase. The underlined terms in the above expression result in a non-analytic momentum dependence of the form p i(µ−λ) , where δ (λ, µ) is a p-independent phase. The the integral can be written as Putting together eq. (4.12) and eq. (4.18), we get up to a constant phase. Going back to the first term in eq. (4.6) with G +− , we observe that the calculation can be repeated with the replacement λ → −λ. However, it can be checked that in this case, the stationary point does not exist for either of the two integrals, and the integrand is highly oscillatory in the entire domain of integration η ∈ (−∞, 0). This gives an exponentially suppressed contribution, and hence can be ignored in our estimate. 6 From this analysis, it becomes clear that the correct way to think of I +− diagrams is fig. 1 (a). In the figure, we represent anti time-ordered part above the line denoting the end of inflation (η ≈ 0) with time running downwards. This is simply a convenient notation to differentiate between I +−/−+ and I ++/−− diagrams. The interaction Hamiltonian evolves the initial Bunch-Davies vacuum into a superposition of excited states |p i ...p n ξ ; q j ...q nχ where n ξ is the number of inflaton excitations while nχ is the number of excitations of the heavy field. The correlation functions essentially picks the amplitude of the state |k 3 ; k 3 from the 'bra' and the amplitude of |k 1 , k 2 ; k 3 from the 'ket'.
The direction of the arrow for chemical potential λ in fig. 1 depends on the phase. For the vertices from time-ordered part, (−η) iλ corresponds to 'injection' of energy into the vertex, while (−η) −iλ corresponds to 'removal' of energy. This is exactly opposite for the anti-time ordered part. This will also help us understand why I ++ /I −− contribution is sub-dominant.
The I ++ diagram is different than I +− in the sense that the produced heavy particle decays into inflaton at a later time as seen in fig. 1 (b). This means that while the energy of order λ is injected at one vertex to boost particle production, similar amount of energy is removed at the other vertex. If λ > µ, the production of massive particle at k 3 vertex is efficient, i.e. stationary point exists for the oscillating integrand. But the decay of the produced particle into inflatons at k 12 vertex is suppressed. This is because the particle has to decay not only into inflatons but also give away energy of order λ which is greater than its energy ∼ µ. This can be seen by the absence of a stationary point for the integral at this vertex, resulting in the suppression of the form e −π(λ−µ) . On the other hand, if µ > λ, Figure 1: Representation of the dominant tree-level contributions from I +− and I ++ diagrams to the bispectrum. For the convenience of notation, the time-ordered part of the diagram is shown below the η ≈ 0 line with time running upwards, while the anti time-ordered part is shown above with time running downwards. An incoming (outgoing) arrow at a vertex denotes injection (extraction) of energy. This makes it plausible that the I +− contribution is enhanced compared to the I ++ contribution, since for the former the "resonant" particle production through energy injection is efficient at both the vertices. Further explanations can be found in the text. the reverse will be true, that is, decay at k 12 will be efficient, but production at k 3 will be suppressed. All in all, we expect I ++ to at least have the suppression of the form e −π|λ−µ| .
We postpone the calculation of the sub-dominant I ++ diagram to section 4.4 (the full analytic treatment can be found in appendix B.2), and instead discuss some of our main results in the following sections.

Constraints
Before looking at the strength of the non-analytic signal, we determine the constraints on various parameters of our model in eq. (3.7). To ensure a controlled derivative expansion in (∂φ) 2 /Λ 4 , we require Λ >φ 1/2 0 . From the Planck data [2], H 4 /φ 2 0 = 8.25 × 10 −8 at k * = 0.05 Mpc −1 for slow-roll inflation. This givesφ 0 ≈ (60H) 2 implying the chemical potential λ 60H. Thus λ can be much larger than the typical energy of the expanding space-time during inflation ∼ H while ensuring EFT control.
We must also ensure that the standard inflaton dynamics is not affected to first order and we have perturbative control over the interactions. The VEV of the complex scalar in eq. (3.14) modifies the kinetic energy of the inflaton as (∂φ This gives a constraint on the VEV as We also require the correction to the power spectrum from the interaction Hamiltonian to be sub-leading. The dominant corrections come from the symmetry-breaking term H int ⊃ αχe −iφ/Λ ∼ −ακ † (ξ/Λ) 2 −iαe −iφ 0 Λ δχ(ξ/Λ)+c.c.. These two contributions are comparable, and can be estimated to be ∆P R /P R ∼ α|κ|/Λ 2 , which should be 1. For λ > M eff , κ ∼ α/λ 2 . Then it is then easy to see than this constraint is stronger than eq. (4.20). Then, the overall constraint is simply α/Λ λ. In order to have only a percent level correction to the power spectrum [2], we take a more conservative value for the coefficient of the symmetry-breaking term, We will use this value of α while evaluating the strength of NG in the next subsection.

Central results
With the technical calculation out of the way, we are now in a position to explore the strength of the non-gaussian signal in different parameter regimes. In this subsection, we present the central results of this paper. Let us estimate the size and the shape of the bispectrum in the squeezed limit using eq. (2.8), where we take Here, we have included the symmetry factor of 2 from k 12 vertex. Inserting values for β 2 , ρ 2 from eqs. (3.30) and (3.32), we get the parametric form of NG in squeezed limit We saturate the constraint on α given in eq. (4.21) and use λ ≡φ 0 /Λ to get the approximate form of NG in the squeezed limit as Let us denote the amplitude of the non-analytic signal in the bispectrum by f oscil such that, where ρ(µ, λ) is some µ and λ-dependent phase. Then, . (4.28) From here onward, we will use |f oscil | to characterize the size of NG.
In fig. 2, we compare above expression from the approximate calculation with the exact calculation from appendix B for a generic value of c = 0.5 and λ = 20, 40. We see that in the region where the approximations are valid, i.e., λ > µ, the parametric form of the non-analytic contribution from the stationary phase analysis is in good agreement with the full analytic form. The small discrepancy can be traced to the contributions from other vertices that we did not take into account, especially the mixing vertex with coefficient β 1 .  Figure 2: Comparison of |f oscil | as a function of M using stationary phase (dashed) and the full analytic form (solid) for two different values of the chemical potential λ = 20H (yellow) and λ = 40H (green), at c = 0.5 for α/Λ = 0.1λ. The parametric dependence matches for µ < λ as expected with a difference of an O(1) factor. This mismatch can be attributed to contributions from other vertices which we did not consider in the stationary phase approximation above just for simplicity. Note the full calculation also smooths out the naive divergence in the stationary phase estimate as µ ≈ λ since for such µ the approximation made in our stationary phase estimate breaks down.
Having confirmed the validity of our approach, we now study the size of NG for various values of the parameters of our model, c and λ, as shown in fig. 3. Some notable features in these plots are: • The most important result is that we obtain an observable NG for masses as large as M ∼ √ cλ. 8 This range is largest for the case of pure chemical potential, i.e, c = 1, as seen in fig. 3(c) where masses as large as λ ( 60H) can be observable.
• The non-analytic component of NG has a weak dependence on M , which gives a plateau-like behavior in |f oscil | for M eff < λ (for c = 0). The signal peaks at M eff λ due to the resonance behaviour. After that, the exponential suppression kicks in and the signal drops as |f oscil | ∼ e −π(µ−λ) .
• For c = 0, M eff is always greater than the chemical potential λ, which means we always suffer the suppression e −π(µ−λ) . Using M eff = √ M 2 + λ 2 for c = 0, we can approximate the suppression for small M as e −π(µ−λ) ≈ e −πM 2 /2λ . M 2 values. In general, the reach for c < 1 can be extended by allowing tachyonic mass for the complex scalar field as long as M 2 eff > 0 is ensured to avoid instability during inflation. For example, the case of c = 0 (with M 2 eff = M 2 + λ 2 ), which appears to have the smallest reach, technically allows O(λ) range of tachyonic masses for the heavy scalar making it effectively equivalent to the case of c = 1. Of course, this equivalence does not hold after the end of inflation where the dynamics would be different for non-tachyonic and tachyonic M 2 whereas the effective mass contribution from c would vanish. It is also interesting to note that the case c = 0 for M 2 > 0 is somewhat analogous to the fermionic chemical potential (see e.g. [26]), where the effective mass is always modified by the chemical potential as M 2 eff = M 2 + λ 2 , since fermions do not have the tachyonic mass or c options. In both cases the reach is limited to M λ 1/2 .
We also demonstrate the oscillatory behavior of the NG signal for a benchmark point M = 5H and λ = 40H in fig. 4. The parameter c affects the frequency of oscillation by modifying M eff as expected. Also, notice that the amplitude of the oscillations has a small dependence on the the ratio k 1 /k 3 , which can be used to infer M eff , as we will see in section 5.
In conclusion, we have demonstrated an enhanced NG signal for heavy complex scalar fields with chemical potential for a wide range of parameter space as seen in fig. 3. The NG is at tree level and hence, the calculation is robust and much simpler than the loop-level processes previously considered in the context of chemical potentials. In fig. 4, we see that the signal has a distinct oscillatory signature with many observable oscillations before the p −3/2 dilution takes over (this dilution is however factored out on the vertical axis in fig. 4).

Sub-dominant contribution: I ++ diagrams
In this section, we go through the calculation of I ++ diagram to convince the reader that it is indeed exponentially suppressed. Taking the same vertices as for the I +− diagram in section 4.1, we write the I ++ contribution There is a time-ordering issue in I ++ /I −− diagrams. But the time scales of the dominant contribution are well separated and thus the integrals can be extended to the entire domain of integration. From the previous calculation of I +− diagrams (we will check for consistency later in the calculation), the integral at k 3 is dominated by earlier time as compared to the integral at k 12 . This simplifies the time ordering. We focus on the part with F ++ propagator which corresponds to fig. 1 (b).
where tilde denotes integrals in I ++ /I −− diagrams.Ĩ (+) k 3 is dominated by the contribution from stationary point at −k 3 η 2 ≡ x * ∼ (λ 2 −µ 2 )/(2λ), and it can be evaluated using similar procedure used to get eq. (4.12), The integral at k 12 vertex on the other hand does not have a stationary point in the domain of integration for λ > µ.
Thenh (x) = − λ x −p+ 1 + µ 2 x 2 is never zero for x > 0 and λ > µ. To evaluate this integral, we instead perform a Wick rotation x → −ix. Then the phase factor e −ipx becomes an exponential e −px which gives contribution in the region x < 1/p. This justifies the timeordering. Since p > 1, we can take the late time limit of the adiabatic mode function which goes as x iµ , Here we have only kept the largest contribution. Taking asymptotic expression for gamma function, Γ(x + iy) y→∞ −−−→ √ 2π|y| x−1/2 e −π|y|/2 , we simplify above expressioñ Combining eq. (4.31) and eq. (4.35), In the regime where λ < µ, we still get the suppression e −π|µ−λ| but from the k 3 vertex instead. We see that I ++ is generally sub-dominant to I +− by a factor of e −π|µ−λ| , and only contributes significantly in a small region when M eff is close to λ.

Symmetry breaking mass term
Till now, we have considered a linear symmetry-breaking term for the complex scalar. However, we can consider a more general derivative coupling of the inflaton to the scalar fields as discussed in eq. (1.5), which can be mapped to a quadratic χχ-type mass correction.
Starting with the form in eq. (1.6), let us write this in terms of a complex scalar χ = σ 1 + iσ 2 , Notice that the first term is a chemical potential ∼ λJ 0 . The second term represents the extent of symmetry-breaking and is ∝ (λ 2 − λ 1 ) as expected. It can be simplified to a symmetry-breaking mass correction using integration by parts, and then rotating the complex field χ → e −iπ/4 χ to get rid of the extra phase factors, In the presence of a Z 2 -symmetry χ → −χ, this constitutes the leading U (1)-breaking. It is also roughly analogous to the studies of chemical potential in fermions, which is another motivation to study its effects in greater detail. The above considerations give rise to a Lagrangian similar to eq. (3.4) but with a quadratic symmetry-breaking m 2 χ 2 instead, Then after the usual field redefinition χ = e −iφ/Λχ , we get We work in the regime H < m M such that the symmetry-breaking does not lead to explosive particle production. The EOM can be written as, Notice that the EOM for χ and χ † are coupled due to the symmetry-breaking mass correction. The diagonalization of non-derivative terms will involve a time-dependent unitary matrix with factors ∝ (−η) ±iλ from the symmetry-breaking mass correction. This will induce off-diagonal elements from the derivative terms in EOM. Invariably, factors of (η) ±iλ can not be removed, which is the source for particle production. Equation (4.42) is similar to narrow parametric resonance with a high frequency mass correction. To get some intuition about the effect of m, we can consider a simplified EOM for a real scalar field σ motivated by eq. (4.42), Given that we will focus on a regime λ > M m > H, we will treat the expansion of spacetime adiabatically compared to particle production. We expect a resonance-like behavior when the frequency of the source ∼ 2λ is twice the natural frequency of the complex scalar, i.e., when ω(η * ) ∼ λ/η * in conformal time. Then the condition for resonance is, (4.44) and it lasts for the duration when the frequency is inside a resonance band. We only consider the first resonance band, where the amplification is the largest [61]: This gives time interval of resonance to be, .  We see that the resonance occurs for a very short duration and hence the amplification is small ∼ O(m 4 /λ 3 ). This is expected since the physical momentum dilutes exponentially fast. Above analysis indicates that if we work in the regime where the amplification is small, i.e., m 4 /λ 3 1, the symmetry-breaking term can be considered as a perturbation. We can then treat it as a Feynman vertex in the "insertion" approximation, while the propagators of the free fields are unchanged. In the following, we estimate its contribution to the bispectrum. The contributions to the bispectrum come from diagrams in fig. 5. The full loop calculation is generally quite challenging, but it can be simplified in certain regimes. In the squeezed limit, the non-analytic signature is prominent when bothχ-modes at the k 3 vertex are soft ∼ k 3 . Then the loop integral can be approximated as With this simplification, the contribution to the bispectrum can be estimated using similar methods as in section 4.1 (except now we have twoχ modes at each vertex and an explicit time-dependence of the form (−η) ±2iλ ), where we have included the loop factor 1/16π 2 . For Λ, µ ∼ O(λ), this expression can be simplified as (4.50) Notice that here the non-analytic signal dilutes faster ∼ p −3 as compared to the case of symmetry-breaking tadpole term, while the non-analytic exponent is twice as large. This is the feature of the contribution being at loop level. The interaction will also modify the power spectrum, which can be expected to be ∼ O(m 4 /λ 3 ). Then to ensure only a percent level correction to the power spectrum, we require (m 4 /λ 3 ) 10 −2 . Using this constraint, |f oscil | can be estimated as, We see that the non-analytic contribution from m 2 χ 2 -type symmetry breaking term is generically somewhat smaller compared to the case of linear symmetry-breaking model as discussed in section 4.3. While the contribution is small, it could still be within the sensitivity of future experiments. Also, this could be the sole contribution if an additional Z 2 -symmetry is imposed on the complex scalar field. One may also consider the parameter regime where the amplification is large and must be included at leading order. We leave a more detailed study of this model to future work.

Inferring effective mass from the bispectrum
In the cosmological collider physics program, the non-analytic exponent of p (= k 12 /k 3 ) in the squeezed limit allows spectroscopy of heavy fields as seen in eq. (1.1). In our case, the dominant contribution to NG has a non-analytic behavior of the form ∼ p −3/2±i(µ−λ) as seen in eq. (4.26), which means an independent measurement of λ is needed to infer the effective mass of the heavy field. In this section, we outline one such procedure which can be used to effectively extract the effective mass and the chemical potential in a certain parameter regime by observing the change in the non-analytic exponent as a function of the squeezing parameter p.
Till now, we have been taking the squeezed limit where p 1, but there is another interesting region when λ p, µ. To get the parametric form of the bispectrum in this regime, we will have to use the following asymptotic form of the hypergeometric function [62], Using this asymptotic form in the analytic expression for the bispectrum (see appendix B for details), we can get the non-analytic dependence on p in the large λ limit, where r is some rational number that depends on the specific interaction. The takeaway is that the asymptotic form in the limit λ p, µ is different than the squeezed limit, and the non-analyticity has the form (1 + p) ±iλ .
This can also be seen from our approximate stationary phase analysis from earlier. For λ > µp, the production at k 12 vertex occurs much earlier when the frequency of the complex scalar is dominated by the physical momentum rather than its mass. Hence, the non-analytic part can not capture the information about M eff , and hence, it is expected to only depend on λ.
In summary, when λ > µp, we expect the non-analytic dependence ∼ (1 + p) ±iλ , and in large p limit, i.e., λ < µp, we expect the non-analyticity of the form ∼ p ±i(µ−λ) . In fig. 6 (a), we see that the full analytic form (blue) agrees with the large λ approximation (yellow) for small values of p. In the squeezed limit on the other hand, the analytic form matches the large p approximation (green). This corroborates our expectation from the previous discussion. (a)  Figure 6: Comparison of the full analytic form evaluated in appendix B (blue) to the large λ approximation (yellow) and large p approximation (green) using (a) the oscillatory form of the dimensionless NG function F and (b) the non-analytic exponent E. We clearly see that as the degree of squeezing is changed the full bispectrum switches between the above two approximations.
To extract the effective mass, we can look at the non-analytic exponent as a function of the degree of squeezing. For this purpose, we define a function E such that E = envelope of (p)∂ p p 3/2 F (k 1 , k 2 , k 3 ) /envelope of p 3/2 F (k 1 , k 2 , k 3 )

(5.4)
Then E is effectively the non-analytic exponent. From fig. 6 (b), we see that it also agrees well with our expectation. Using this feature, we can extract µ = M 2 eff − 9/4 (and hence M eff ) as well as the chemical potential λ. This procedure is more efficient for c close to 1, i.e, the case of pure chemical potential. It becomes harder to extract the effective mass if c 1, when µ ∼ O(λ), and we are effectively always in the large p limit.

Conclusions
In this paper, we have shown that a derivative coupling of the inflaton to heavy scalar fields of "chemical potential" type can lead to unsuppressed production of heavy particles and leave observable imprints in primordial non-Gaussianities. In particular, the chemical potential can overturn the Boltzmann suppression for masses far above the inflationary Hubble scale H, ordinarily present in cosmological particle production. This effect is quite general, requiring only two (or more) real scalar fields where the symmetry rotating one into the other is broken. In essence, the unsuppressed particle production above H draws on the large kinetic energy of the inflaton background,φ 0 ∼ (60H) 2 . A heavy mass reach up to ∼ 60H can thereby be accomplished while still remaining within a controlled derivative expansion for the inflationary effective field theory. Signatures of such heavy masses could be observed in the bispectrum with strength f NL ∼ O(0.01 − 10), within the sensitivity of upcoming LSS and future 21-cm experiments.
The main results of our analysis are summarized in fig. 3 and fig. 4. In our model, the contribution to the bispectrum is at tree-level, which removes the loop-suppression seen in earlier spin-1/2 and spin-1 chemical potential examples. Importantly, it also makes the calculations analytically tractable and physically transparent. We outlined a procedure to infer the effective mass of the heavy field, which relies on the variation in the non-analytic exponent of the co-moving momentum ratio (degree of "squeezing") as in fig. 6.
It is straightforward to generalize our mechanism to more than one heavy complex scalar, each with its own distinct chemical potential λ, thereby allowing us to explore multiple "channels". More non-trivially, our mechanism can be extended to massive complex fields with arbitrary non-zero spins. The mechanism can also be extended to the curvaton paradigm [63][64][65] for inflation, which can exhibit new features. We will take up these generalizations in future work.
where f (x) = e πµ/2 x 3/2 H (2) iµ (x) is a solution to the EOM for mode functions as seen in section 3.4. Then, From here, we can systematically move the integration past the differential operator using integration by parts to get a differential equation for I 2 with respect to variable p. To clarify the point, consider the last term, Similarly for the derivative terms, we integrate by parts. For example, the second term can be written as, The first term is the boundary term which can be dropped (justify). Then we are left with, Following this procedure, we get a differential equation for I 2 (n, p) as, (p 2 − 1)∂ 2 p + 2p n + 3 2 ∂ p + M 2 + n − 1 2 n + 5 2 I 2 (n, p) = 0 . (A.7) We want the solution to be regular at p = 1. This implies, I 2 (n, p) ∝ 2 F 1 n + 1 + iµ, n + 1 − iµ, n + 3/2; (1 − p) 2 . (A.8) The normalization is fixed by explicitly carrying out the integration in large p limit using small x expansion of the Hankel function, and then matching it to the large p expansion of hypergeometric function. We also note that this derivation holds even if n has an imaginary part.