BINGO: A code for the efficient computation of the scalar bi-spectrum

# BINGO: A code for the efficient computation of the scalar bi-spectrum

Dhiraj Kumar Hazra Harish-Chandra Research Institute, Chhatnag Road, Jhunsi, Allahabad 211019, India111Current affiliation: Asia Pacific Center for Theoretical Physics, Pohang, Gyeongbuk 790-784, Korea. E-mail: dhiraj@apctp.org..    L. Sriramkumar Department of Physics, Indian Institute of Technology Madras, Chennai 600036, India.    Jérôme Martin Institut d’Astrophysique de Paris, UMR7095-CNRS, Université Pierre et Marie Curie, 98bis boulevard Arago, 75014 Paris, France.
July 14, 2019
###### Abstract

We present a new and accurate Fortran code, the BI-spectra and Non-Gaussianity Operator (BINGO), for the efficient numerical computation of the scalar bi-spectrum and the non-Gaussianity parameter in single field inflationary models involving the canonical scalar field. The code can calculate all the different contributions to the bi-spectrum and the parameter for an arbitrary triangular configuration of the wavevectors. Focusing firstly on the equilateral limit, we illustrate the accuracy of BINGO by comparing the results from the code with the spectral dependence of the bi-spectrum expected in power law inflation. Then, considering an arbitrary triangular configuration, we contrast the numerical results with the analytical expression available in the slow roll limit, for, say, the case of the conventional quadratic potential. Considering a non-trivial scenario involving deviations from slow roll, we compare the results from the code with the analytical results that have recently been obtained in the case of the Starobinsky model in the equilateral limit. As an immediate application, we utilize BINGO to examine of the power of the non-Gaussianity parameter to discriminate between various inflationary models that admit departures from slow roll and lead to similar features in the scalar power spectrum. We close with a summary and discussion on the implications of the results we obtain.

###### pacs:
98.80.Cq, 98.70.Vc

## 1 Non-Gaussianities and primordial features

Over the last half-a-dozen years, it has been increasingly realized that the detection of non-Gaussianities in the primordial perturbations can considerably help in constraining the inflationary models (see Refs. [1, 2, 3]; for early efforts in this direction, see Refs. [4]). In particular, the detection of a high value for the dimensionless non-Gaussianity parameter that is often used to describe the amplitude of the scalar bi-spectrum can rule out a wide class of models. For instance, if the extent of non-Gaussianity actually proves to be as large as the mean values of arrived at from the WMAP data (see Refs. [5, 6, 7]; in this context, also see Refs. [8, 9]), then canonical scalar field models that lead to slow roll inflation and nearly scale invariant primordial spectra will cease to be consistent with the data.

However, the evaluation of the scalar bi-spectrum in a generic inflationary model proves to be a non-trivial task [10]. Often, it is the slow roll approximation that is resorted to in order to arrive at analytical expressions for the scalar bi-spectrum and the non-Gaussianity parameter  [1, 2]. While the slow roll approximation can actually encompass a relatively wide class of models, evidently, the approach will cease to be valid when departures from slow roll occur. In such situations, one is often forced to approach the problem numerically. In this work, we shall present a new and accurate Fortran code, the BI-spectra and Non-Gaussianity Operator (BINGO), to numerically evaluate the scalar bi-spectrum and the non-Gaussianity parameter  for single field inflationary models involving the canonical scalar field. Although, some partial numerical results have already been published in the literature, we believe that it is for the first time that a general and efficient code has been put together to evaluate the complete scalar bi-spectrum. Though, we shall largely restrict our discussion in this paper to the equilateral case, as we shall illustrate, the code can compute the scalar bi-spectrum and the parameter  for any triangular configuration of the wavevectors. Also, as we shall demonstrate, BINGO can also compute all the different contributions to the bi-spectrum. Moreover, it is worth noting that, under certain conditions, the code can arrive at the results within a matter of a few minutes. We should mention here that we have made a limited version of BINGO, viz. one that focuses on the equilateral limit, available online at https://www.physics.iitm.ac.in/~sriram/bingo/bingo.html.

As an immediate application, we shall utilize the code to examine of the power of the non-Gaussianity parameter to discriminate between different inflationary models that admit deviations from slow roll and generate similar features in the scalar power spectrum. Recall that, most single field inflationary models naturally lead to an extended period of slow roll and hence to nearly scale invariant primordial power spectra (see, for example, any of the following texts [11] or reviews [12]), which seem to be fairly consistent with the recent data on the anisotropies in the Cosmic Microwave Background (CMB) (in this context, see Refs. [5, 6, 7]; for a comparison of a class of inflationary models with the data, see, for example, Ref. [13] and references therein) as well as other observational constraints. However, it has been repeatedly noticed that certain features in the inflationary scalar power spectrum can improve the fit to the CMB data at the cost of some additional parameters (see, for instance, Refs. [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28]). Though the statistical significance of these features remain to be understood satisfactorily [29], they gain importance from the phenomenological perspective of comparing the models with the data, because only a smaller class of single field inflationary models, which allow for departures from slow roll, can generate them. Interestingly, demanding the presence of features in the scalar power spectrum seems to generically lead to larger non-Gaussianities (see, for example, Refs. [30, 31]). Therefore, features may offer the only route (unless one works with non-vacuum initial states [32]) for the canonical scalar fields to remain viable if turns out to be significant.

If indeed the presence of features turns out to be the correct reason behind possibly large non-Gaussianities, can we observationally identify the correct underlying inflationary scenario, in particular, given the fact that different models can lead to similar features in the scalar power spectrum? In other words, to what extent can the non-Gaussianity parameter help us discriminate between the inflationary models that permit features? To address this question, we shall consider a few typical inflationary models leading to features, assuming that they can be viewed as representatives of such a class of scenarios. Concretely, we shall consider the Starobinsky model [33] and the punctuated inflationary scenario [17], both of which result in a sharp drop in power at large scales that is followed by oscillations. We shall also study large and small field models with an additional step introduced in the inflaton potential [19, 30, 31]. The step leads to a burst of oscillations in the scalar power spectrum which improve the fit to the outliers near the multipole moments of and in the CMB anisotropies. We shall also consider oscillating inflaton potentials such as the one that arises in the axion monodromy model which lead to modulations in the power spectrum over a wide range of scales [25, 26, 27, 28].

The plan of this paper is as follows. In the following section, we shall quickly describe a few essential details pertaining to the power spectrum and the bi-spectrum as well as the non-Gaussianity parameter in inflationary models involving a single, canonical, scalar field. In Sec. 3, after demonstrating that the super Hubble contributions to the bi-spectrum during inflation prove to be negligible, we shall describe the method that BINGO adopts to numerically compute the bi-spectrum and the non-Gaussianity parameter . We shall also illustrate the extent of accuracy of BINGO by comparing the numerical results with the analytical results available in three situations. Firstly, restricting ourselves to the equilateral limit, we shall compare the results from the numerical computations with the spectral dependence of the bi-spectrum expected in power law inflation. Secondly, focusing on the case of the archetypical quadratic potential, we shall contrast the numerical results with the analytical expressions obtained using the slow roll approximation for an arbitrary triangular configuration of the wavevectors. Thirdly, considering a more non-trivial situation involving a brief period of fast roll, we shall compare the numerical results with the analytical results that have recently been obtained in the case of the Starobinsky model in the equilateral limit [34, 35]. In the succeeding two sections, we shall utilize BINGO to examine of the power of the non-Gaussianity parameter to discriminate between different inflationary models that permit deviations from slow roll and generate similar features in the scalar power spectrum. After a quick outline of the inflationary models of our interest, we shall discuss the scalar power spectra that arise in these models. We shall then present the main results, and compare the that arise in the various models that we consider. We shall finally conclude in Sec. 6 with a brief summary and outlook.

A few words on our conventions and notations are in order before we proceed. We shall work with units such that , and we shall set . As is often done in the context of inflation, we shall assume the background to be described by the spatially flat, Friedmann-Lemaitre-Robertson-Walker line-element. We shall denote the conformal time as , while shall represent the number of e-folds. Moreover, an overprime shall denote differentiation with respect to the conformal time coordinate . Lastly, shall denote the scale factor, with the corresponding Hubble parameter being defined as .

## 2 The scalar power spectrum, bi-spectrum and the parameter fNL

In this section, we shall first rapidly summarize the essential definitions of the power spectrum and the bi-spectrum as well as the non-Gaussianity parameter . We shall then list the various contributions to the bi-spectrum that arise in the Maldacena formalism (for the original discussions, see Refs. [1, 2, 3, 30, 31]; for more recent discussions, see, for example, Refs. [34, 36]).

### 2.1 Essential definitions and relations

Let us begin by recalling some essential points concerning the power spectrum. On quantization, the operator corresponding to the curvature perturbation can be expressed as

 ^R(η,x) = ∫d3k(2π)3/2^Rk(η)eik⋅x (1) = ∫d3k(2π)3/2[^akfk(η)eik⋅x+^a†kf∗k(η)e−ik⋅x],

where and are the usual creation and annihilation operators that satisfy the standard commutation relations. The modes are governed by the differential equation [11, 12]

 f′′k+2z′zf′k+k2fk=0, (2)

where , with being the first slow roll parameter. The dimensionless scalar power spectrum is defined in terms of the correlation function of the Fourier modes of the curvature perturbation as follows:

 ⟨0|^Rk(η)^Rp(η)|0⟩=(2π)22k3PS(k)δ(3)(k+p), (3)

where denotes the Bunch-Davies vacuum which is defined as  [37]. In terms of the modes , the scalar power spectrum is given by

 PS(k)=k32π2|fk|2. (4)

As is well known, numerically, the initial conditions are imposed on the modes when they are well inside the Hubble radius, and the power spectrum is evaluated at suitably late times when the modes are sufficiently outside the Hubble radius (see, for instance, Refs. [38]). We shall discuss the details concerning the numerical evolution of the background as well as the perturbations and the computation of the corresponding power spectrum in the next section.

The scalar bi-spectrum is related to the three point correlation function of the Fourier modes of the curvature perturbation , evaluated at the end of inflation, say, at the conformal time , as follows [6]:

 ⟨^Rk1(ηe)^Rk2(ηe)^Rk3(ηe)⟩=(2π)3BS(k1,k2,k3)δ(3)(k1+k2+k3). (5)

In our discussion below, for the sake of convenience, we shall set

 BS(k1,k2,k3)=(2π)−9/2G(k1,k2,k3). (6)

The observationally relevant non-Gaussianity parameter  is introduced through the relation

 (7)

where denotes the Gaussian quantity, and the factor of arises due to the relation between the Bardeen potential and the curvature perturbation during the matter dominated epoch. Upon making use of the corresponding relation between and in Fourier space and the Wick’s theorem, one obtains that [1, 2, 3]

 ⟨^Rk1^Rk2^Rk3⟩ = −3fNL10(2π)4(2π)−3/21k31k32k33δ(3)(k1+k2+k3) (8) ×[k31PS(k2)PS(k3)+two permutations].

This expression can then be utilized to arrive at the following relation between the non-Gaussianity parameter and the scalar bi-spectrum or, equivalently, the quantity that we have introduced [34, 36]:

 fNL(k1,k2,k3) = −103(2π)−4(k1k2k3)3G(k1,k2,k3) (9) ×[k31PS(k2)PS(k3)+two permutations]−1.

In particular, in the equilateral limit, i.e. when , this expression simplifies to

 feqNL(k)=−1091(2π)4k6G(k)P2S(k). (10)

Two points need to be emphasized here regarding the above expressions for the quantities and . Firstly, it ought to be noted that the non-Gaussianity parameter introduced in Eq. (7) [as well as the expression arrived at in Eq. (9)] corresponds to the so-called local form of the parameter. Secondly, the quantity corresponds to the equilateral limit of the local parameter, and is intrinsically different from a similar parameter introduced when the bi-spectrum has the equilateral form (in this context, see, for instance, Ref. [7]).

### 2.2 The scalar bi-spectrum in the Maldacena formalism

In the Maldacena formalism [1], the bi-spectrum is evaluated using the standard rules of perturbative quantum field theory, based on the interaction Hamiltonian that depends cubically on the curvature perturbation. It can be shown that the bi-spectrum that results from the interaction Hamiltonian can be expressed as [1, 2, 3, 30, 31, 34, 36]

 G(k1,k2,k3) ≡ 7∑C=1GC(k1,k2,k3) (11) ≡ +G7(k1,k2,k3).

The quantities with correspond to the six terms in the interaction Hamiltonian, and are described by the integrals

 G1(k1,k2,k3) = 2i∫ηeηidηa2ϵ21(f∗k1f′∗k2f′∗k3+two permutations), (12) G2(k1,k2,k3) = −2i(k1⋅k2+two permutations) (13) ×∫ηeηidηa2ϵ21f∗k1f∗k2f∗k3, G3(k1,k2,k3) = −2i∫ηeηidηa2ϵ21[(k1⋅k2k22)f∗k1f′∗k2f′∗k3 (14) +five permutations], G4(k1,k2,k3) = i∫ηeηidηa2ϵ1ϵ′2(f∗k1f∗k2f′∗k3+two permutations), (15) G5(k1,k2,k3) = i2∫ηeηidηa2ϵ31[(k1⋅k2k22)f∗k1f′∗k2f′∗k3 (16) +five permutations], G6(k1,k2,k3) = i2∫ηeηidηa2ϵ31{[k21(k2⋅k3)k22k23]f∗k1f′∗k2f′∗k3 (17) +two permutations},

where is the second slow roll parameter that is defined with respect to the first as follows: . The lower limit of the above integrals, viz. , denotes a sufficiently early time when the initial conditions are imposed on the modes. The additional, seventh term arises due to the field redefinition, and its contribution to is given by

 G7(k1,k2,k3)=ϵ2(ηe)2(|fk2(ηe)|2|fk3(ηe)|2+two permutations). (18)

## 3 The numerical computation of the scalar bi-spectrum

In this section, after illustrating that the super-Hubble contributions to the complete bi-spectrum during inflation proves to be negligible, we shall outline the methods that BINGO adopts to numerically evolve the equations governing the background and the perturbations, and eventually evaluate the inflationary scalar power and bi-spectra. Also, we shall illustrate the extent of accuracy of BINGO by comparing the results from the code with the expected form of the bi-spectrum in the equilateral limit in power law inflation and the analytical results that are available in the case of the Starobinsky model [34, 35]. We shall also contrast the numerical results with the analytical results obtained under the slow roll approximation in the case of the popular quadratic potential for an arbitrary triangular configuration of the wavevectors.

### 3.1 The contributions to the bi-spectrum on super-Hubble scales

It is clear from the expressions in the previous section that the evaluation of the bi-spectrum involves integrals over the mode and its derivative as well as the slow roll parameters , and the derivative . While evaluating the power spectra, it is well known that, in single field inflationary models, it suffices to evolve the curvature perturbation from an initial time when the modes are sufficiently inside the Hubble radius to a suitably late time when the amplitude of the curvature perturbation settles down to a constant value on super-Hubble scales (see, for example, Refs. [38]). We shall illustrate that many of the contributions to the bi-spectrum prove to be negligible when the modes evolve on super-Hubble scales. Interestingly, we shall also show that, those contributions to the bi-spectrum which turn out to be significant at late times when the modes are well outside the Hubble radius are canceled by certain other contributions that arise. As a consequence, we shall argue that, numerically, it suffices to evaluate the integrals over the period of time during which the curvature perturbations have been conventionally evolved to arrive at the power spectra, viz. from the sub-Hubble to the super-Hubble scales. In fact, we had recently shown that, for cosmological modes that are on super-Hubble scales during preheating, the contributions to the bi-spectrum and, equivalently, to the non-Gaussianity parameter , prove to be negligible [36]. The calculations we had presented in the work can be easily extended to the case of the modes on super-Hubble scales during inflation. Therefore, we shall rapidly sketch the essentials arguments below, referring the readers to our recent work for further details [36].

#### 3.1.1 Evolution of fk on super-Hubble scales

During inflation, when the modes are on super-Hubble scales, it is well known that the solution to can be written as [11, 12]

 fk≃Ak+Bk∫ηd~ηz2(~η), (19)

where and are -dependent constants which are determined by the initial conditions imposed on the modes in the sub-Hubble limit. The first term involving is the growing mode, which is actually a constant, while the term containing represents the decaying mode. Therefore, on super-Hubble scales, the mode simplifies to

 fk≃Ak. (20)

Moreover, the leading non-zero contribution to its derivative is determined by the decaying mode, and is given by

 f′k≃Bkz2=¯Bka2ϵ1, (21)

where we have set .

It is now a matter of making use of the above solutions for and to determine the super-Hubble contributions to the bi-spectrum during inflation.

#### 3.1.2 The various contributions

To begin with, note that, each of the integrals , where , can be divided into two parts as follows:

 GC(k1,k2,k3)=GisC(k1,k2,k3)+GseC(k1,k2,k3). (22)

The integrals in the first term run from the earliest time (i.e. ) when the smallest of the three wavenumbers , and is sufficiently inside the Hubble radius [typically corresponding to ] to the time (say, ) when the largest of the three wavenumbers is well outside the Hubble radius [say, when ]. Then, evidently, the second term will involve integrals which run from the latter time to the end of inflation at  (in this context, see Fig. 1).

In what follows, we shall discuss the various contributions to the bi-spectrum due to the terms . We shall show that the corresponding contribution either remains small or, when it proves to be large, it is exactly canceled by another contribution to the bi-spectrum.

Let us first focus on the fourth term since it has often been found to lead to the largest contribution to the bi-spectrum when deviations from slow roll occur [30, 31, 34, 36]. As the slow roll parameters turn large towards the end of inflation, we can expect this term to contribute significantly at late times. However, as we shall quickly illustrate, such a late time contribution is exactly canceled by the contribution from which arises due to the field redefinition [36]. Upon using the form (20) of the mode and its derivative (21) on super-Hubble scales in the expression (15), the integral involved can be trivially carried out with the result that the corresponding contribution to the bi-spectrum can be expressed as

 Gse4(k1,k2,k3) ≃ iM2Pl[ϵ2(ηe)−ϵ2(ηs)][|Ak1|2|Ak2|2(Ak3¯B∗k3−A∗k3¯Bk3) (23) +two permutations].

The Wronskian corresponding to the equation of motion (2) and the standard Bunch-Davies initial condition leads to the relation , which can then be utilized to arrive at the following simpler expression [36]:

 Gse4(k1,k2,k3) ≃ −12[ϵ2(ηe)−ϵ2(ηs)] (24) ×[|Ak1|2|Ak2|2+ two permutations].

The first of these terms involving the value of at the end of inflation exactly cancels the contribution [with set to in Eq. (18)] that arises due to the field redefinition. But, the remaining contribution cannot be ignored and needs to be taken into account. It is useful to note that this term is essentially the same as the one due to the field redefinition, but which is now evaluated on super-Hubble scales (i.e. at ) rather than at the end of inflation. In other words, if we consider the fourth and the seventh terms together, it is equivalent to evaluating the contribution to the bi-spectrum corresponding to , and adding to it the contribution due to the seventh term evaluated at rather than at the end of inflation.

Let us now turn to the contribution due to the second term, which can occasionally prove to be comparable to the contribution due to the fourth term [34]. Upon making use of the behavior of the mode and its derivative on super-Hubble scales in the integral (13), one obtains the contribution to the bi-spectrum due to to be [34, 36]

 Gse2(k1,k2,k3) = −2iM2Pl(k1⋅k2+two permutations) (25) ×|Ak1|2|Ak2|2|Ak3|2[I2(ηe,ηs)−I∗2(ηe,ηs)],

where the quantity is described by the integral

 I2(ηe,ηs)=∫ηeηsdηa2ϵ21. (26)

Note that, due to the quadratic dependence on the scale factor, actually, is a rapidly growing quantity at late times. However, the complete super-Hubble contribution to the bi-spectrum vanishes identically since the integral is a purely real quantity [36]. Hence, in the case of the second term, it is sufficient to evaluate the contribution to the bi-spectrum due to .

Due to their structure, one finds that the first and the third terms and the fifth and the sixth terms can be evaluated together. On super-Hubble scales, one can easily show that the contributions due to the first and the third terms can be written as

 Ges1(k1,k2,k3)+Ges3(k1,k2,k3) = 2iM2Pl[(1−k1⋅k2k22−k1⋅k3k23)|Ak1|2 (27) ×(Ak2¯B∗k2Ak3¯B∗k3−A∗k2¯Bk2A∗k3¯Bk3) + two permutations]I13(ηe,ηs),

while the corresponding contributions due to the fifth and the sixth terms can be obtained to be

 Gse5(k1,k2,k3)+Gse6(k1,k2,k3) = iM2Pl2 (28) ×[(k1⋅k2k22+k1⋅k3k23+k21(k2⋅k3)k22k23) ×|Ak1|2 ×(Ak2¯B∗k2Ak3¯B∗k3−A∗k2¯Bk2A∗k3¯Bk3) + two permutations]I56(ηe,ηs),

where the quantities and are described by the integrals

 I13(ηe,ηs)=∫ηeηsdηa2 (29)

and

 I56(ηe,ηs)=∫ηeηsdηa2ϵ1. (30)

Hence, the non-zero, super-Hubble contribution to the bi-spectrum is determined by the complete contribution due to the first, the third, the fifth and the sixth terms arrived at above. In order to illustrate that this contribution is insignificant, we shall now turn to estimating the amplitude of the corresponding contribution to the non-Gaussianity parameter .

### 3.2 An estimate of the super-Hubble contribution to the non-Gaussianity parameter

Let us restrict ourselves to the equilateral limit for simplicity. In such a case, the super-Hubble contributions due to the first, the third, the fifth and the sixth terms to the non-Gaussianity parameter  can be obtained by substituting the quantities (27) and (28) above in the expression (10). It is straightforward to show that the corresponding is given by

 feq(se)NL(k) ≃ −5iM2Pl18⎛⎝A2k¯B∗k2−A∗k2¯B2k|Ak|2⎞⎠ (31) ×[12I13(ηe,ηs)−94I56(ηe,ηs)],

where we have made use of the fact that at late times in order to arrive at the power spectrum.

To estimate the above super-Hubble contribution to the non-Gaussianity parameter , let us focus on inflation of the power law form. In power law inflation, the scale factor is given by

 a(η)=a1(ηη1)γ+1, (32)

with and being constants, while is a free index. It is useful to note that, in such a case, the first slow roll parameter is a constant, and is given by . Under these conditions, the quantity within the brackets involving and in the expression above can be easily evaluated (in this context, see Ref. [36]) to arrive at

 feq(se)NL(k) = 572π[12−9(γ+2)γ+1]Γ2(γ+12)22γ+1(2γ+1)(γ+2) (33) ×(γ+1)−2(γ+1)sin(2πγ)[1−HsHee−3(Ne−Ns)] ×(kasHs)−(2γ+1).

It should be mentioned that, in arriving at this expression, for convenience, we have set to be , which corresponds to being , viz. the scale factor at . Moreover, while and denote the e-folds corresponding to and , and represent the Hubble scales at these times, respectively. Recall that, denotes the conformal time when the largest wavenumber of interest, say, , is well outside the Hubble radius, i.e. when . Since is expected to be at least for the smallest cosmological scale, it is clear that the factor involving can be completely neglected. Observations point to the fact that . Therefore, if we further assume that , where , we find that the above estimate for the non-Gaussianity parameter reduces to

 feq(se)NL(k)≲−5ε29(ksasHs)3≃−10−19, (34)

where, in obtaining the final value, we have set . The inequality above arises due to the fact that, for larger scales, i.e. when , at . In models involving the canonical scalar field, the smallest values of are typically generated in slow roll inflationary scenarios, wherein the non-Gaussianity parameter has been calculated to be of the order of the first slow roll parameter [1, 3]. The above estimate clearly points to fact that the super-Hubble contributions to the complete bi-spectrum and the non-Gaussianity parameter can be entirely ignored.

In summary, to determine the scalar bi-spectrum, it suffices to evaluate the contributions to the bi-spectrum due to the quantities , with , which involve integrals running from the initial time to the time when the smallest of the three modes reaches super-Hubble scales (cf. Fig. 1). Further, the addition of the contribution due to the field redefinition evaluated at ensures that no non-trivial super-Hubble contributions are ignored. In the following sub-section, with the help of a specific example, we shall also corroborate these conclusions numerically.

### 3.3 Details of the numerical methods

The scalar bi-spectrum and the parameter can be easily evaluated analytically in the slow roll inflationary scenario [1]. However, barring some exceptional cases [34, 35, 39], it often proves to be difficult to evaluate the bi-spectrum analytically when departures from slow roll occur. Hence, one has to resort to numerical computations in such cases.

BINGO solves the background as well as the perturbation equations using RKSUITE [40], which is a publicly available routine to solve ordinary differential equations. We shall treat the number of e-folds as the independent variable, which allows for an efficient and accurate computation. To obtain the power spectrum, we impose the standard Bunch-Davies initial conditions on the perturbations when the modes are well inside the Hubble radius, and evolve them until suitably late times. Typically, in the case of smooth inflaton potentials, it suffices to evolve the modes from an initial time when . However, in certain models wherein the potentials contain oscillatory terms, the modes may have to be evolved from deeper inside the Hubble radius. For example, in the case of the axion monodromy model that we shall discuss in due course (see Sub-sec. 4.2), for values of the model parameters of our interest, we find that we need to evolve the modes from an initial time when , so that the resonance that occurs in the model due to the oscillations in the potential is captured [26, 28, 30, 31]. The scalar power spectra that arise in the various models of our interest, as displayed in Fig. 9, have all been evaluated at super-Hubble scales, say, when , which is typically when the amplitude of the curvature perturbations freeze in.

Having obtained the behavior of the background and the modes, BINGO carries out the integrals involved in arriving at the bi-spectrum using the method of adaptive quadrature [41]. It is useful to note that, in the equilateral limit of the bi-spectrum, which we shall largely concentrate on in this work, we can evolve each of the modes of interest independently and calculate the integrals for the modes separately. The integrals actually contain a cut off in the sub-Hubble limit, which is essential for singling out the perturbative vacuum [1, 2, 3]. Numerically, the presence of the cut off is fortunate since it controls the contributions due to the continuing oscillations that would otherwise occur. We should highlight here that such a cut off procedure was originally introduced in the earliest numerical efforts to evaluate the scalar bi-spectrum in situations involving deviations from slow roll (in this context, see Refs. [30]). Generalizing the cut off that is often introduced analytically in the slow roll case, we shall impose a cut off of the form , where is a small parameter. In the previous two sub-sections, we had discussed as to how the integrals need to be carried out from the early time when the largest scale is well inside the Hubble radius to the late time when the smallest scale is sufficiently outside. If one is focusing on the equilateral configuration, rather than integrate from to , it suffices to compute the integrals for the modes from the time when each of them satisfy the sub-Hubble condition, say, , to the time when they are well outside the Hubble radius, say, when . In other words, one carries out the integrals exactly over the period the modes are evolved to obtain the power spectrum. The presence of the cut-off ensures that the contributions at early times, i.e. near , are negligible. Furthermore, it should be noted that, in such a case, the corresponding super-Hubble contribution to will saturate the bound (34) in power law inflation for all the modes.

With the help of a specific example, let us now illustrate that, for a judicious choice of , the results from BINGO are largely independent of the upper and the lower limits of the integrals. In fact, working in the equilateral limit, we shall demonstrate these points in two steps for the case of the standard quadratic potential [see Eq. (46)]. Firstly, focusing on a specific mode, we shall fix the upper limit of the integral to be the time when . Evolving the mode from different initial times, we shall evaluate the integrals involved from these initial times to the fixed final time for different values of . This exercise helps us to identify an optimal value for when we shall eventually carry out the integrals from . Secondly, upon choosing the optimal value for and integrating from , we shall calculate the integrals for different upper limits. For reasons outlined in the previous two sub-sections, it proves to be necessary to consider the contributions to the bi-spectrum due to the fourth and the seventh terms together. Moreover, since the first and the third, and the fifth and the sixth, have similar structure in the equilateral limit, it turns out to be convenient to club these terms as we have discussed before. In Fig. 2, we have plotted the value of times the different contributions to the bi-spectrum, viz. , , and , as a function of when the integrals have been carried out from of , and for a mode which leaves the Hubble radius around e-folds before the end of inflation in the case of the quadratic potential. The figure suggests that the term is fairly independent of the cut-off parameter . This can be attributed to the fact that the integral depends on the quantity , which effectively behaves as a cut off. For this reason, actually, we shall not introduce the cut off while evaluating the term . Moreover, the figure clearly indicates to be a highly suitable value for the other terms. A larger leads to a sharper cut off reducing the value of the integrals. One could work with a smaller , in which case, the figure suggests that, one would also need to necessarily integrate from deeper inside the Hubble radius.

In Fig. 3, after fixing to be and, with the initial conditions imposed at , we have plotted the four contributions to the bi-spectrum for a mode that leaves the Hubble radius at e-folds before the end of inflation as a function of the upper limit of the integrals.

It is evident from the figure that the values of the integrals converge quickly once the mode leave the Hubble radius. For efficient numerical integration, as in the case of the power spectrum, we have chosen the super-Hubble limit to correspond to . We have repeated similar tests for the other models of our interest too. These tests confirm the conclusions that we have arrived at above, indicating the robustness of the numerical methods and procedures that we have adopted. It is important that we stress three points further at this stage of our discussion. Firstly, in situations involving departures from slow roll, the dominant contribution, viz. the term, does not require the introduction of the cut off at all. Secondly, in the absence of analytical expressions to compare the numerical results with, it is imperative that a suitable value for the cut off parameter is arrived at by repeating the exercise that we have carried out in Fig. 2 for the specific model of interest. In this context, it is further important to appreciate the fact that the integrals involved will depend on both the cut off parameter as well as the initial value of . A complete check seems to be mandatory in order to ensure that a unique result is arrived at. (As we mentioned above, we have indeed carried out this exercise for the models that we consider here.) Lastly, we also need to emphasize that the above comments and conclusions would not apply to certain special cases wherein the amplitude of the curvature perturbation can evolve on super-Hubble scales such as, say, in ultra slow roll inflation (in this context, see, for instance, Ref. [42]).

Before we go on to consider the bi-spectra generated in the inflationary models of our interest, in order to demonstrate the accuracy of BINGO, we shall compare the numerical results from the code with the analytical results that can be arrived at in three situations. In the following sub-sections, we shall compare the results from BINGO with: (i) the spectral dependence that can be arrived at in the context of power law inflation in the equilateral limit, (ii) the slow roll results as applied to the case of the conventional quadratic potential for a generic triangular configuration of the wavevectors, and (iii) the analytic results that can be obtained for the Starobinsky model in the equilateral limit.

### 3.4 Comparison with the analytical results in the case of power law inflation

The first case that we shall discuss is power law inflation wherein, as we shall soon outline, the spectral shape of the non-zero contributions to the bi-spectrum can be easily arrived at in the equilateral limit.

Consider the case of power law inflation described by the scale factor (32) with . In such a case, as we have seen, is a constant and, hence, and , which involve derivatives of , reduce to zero. Since the contributions due to the fourth and the seventh terms, viz. and , depend on and , respectively [cf. Eqs. (15) and (18)], these terms vanish identically in power law inflation. As is well known, one can express the modes as , where is the Mukhanov-Sasaki variable. In power law inflation, one finds that, the variable depends only on the combination (see, for example, Ref. [36]). Moreover, since is a constant in power law expansion, we have . Under these conditions, in the equilateral limit, with a simple rescaling of the variable of integration in the expressions (12), (13), (14), (16) and (17), it is straightforward to show that the quantities , , , and , all depend on the wavenumber as . Then, upon making use of the asymptotic form of the modes , it is easy to illustrate that the corresponding contributions to the bi-spectrum, viz. , and , all behave as . Since the power spectrum in power law inflation is known to have the form (see, for example, Refs. [43]), the expression (10) for then immediately suggests that the quantity will be strictly scale invariant for all . In fact, apart from these results, it is also simple to establish the following relation between the different contributions: , a result, which, in fact, also holds in slow roll inflation [34]. In other words, in power law inflation, it is possible to arrive at the spectral dependence of the non-zero contributions to the bi-spectrum without having to explicitly calculate the integrals involved. Further, one can establish that the non-Gaussianity parameter is exactly scale independent for any value of . While these arguments do not help us in determining the amplitude of the various contributions to the bi-spectrum or the non-Gaussianity parameter, their spectral shape and the relative magnitude of the above-mentioned terms provide crucial analytical results to crosscheck our numerical code. In Fig. 4, we have plotted the different non-zero contributions to the bi-spectrum computed using our numerical code and the spectral dependence we have arrived at above analytically for two different values of  in the case of power law inflation. We have also indicated the relative magnitude of the first and the third and the fifth and the sixth terms arrived at numerically. Lastly, we have also illustrated the scale independent behavior of the non-Gaussianity parameter for both the values of . It is clear from the figure that the numerical results agree well with the results and conclusions that we arrived at above analytically.

### 3.5 Comparison for an arbitrary triangular configuration

Let us now turn to the example of the conventional quadratic potential. The bi-spectrum in such a case can be evaluated analytically in the slow roll approximation for an arbitrary triangular configuration of the wavevectors [1, 2, 3]. In the slow roll limit, one finds that, the contributions due to the first, the second, the third and the seventh terms are of the same order, while the remaining terms turn out to be comparatively insignificant. For convenience, let us explicitly write down the contributions to the bi-spectrum due to the first, the second, the third and the seventh terms arrived at in the slow roll approximation. They are given by

 G1(k1,k2,k3) = H4I16M4Plϵ1(1k1k2k3)3 (35) ×[(1+k1kT)k22k23kT+two permutations], G2(k1,k2,k3) = H4I16M4Plϵ1(1k1k2k3)3(k1⋅k2+two permutations) (36) ×[−kT+1kT(k