{'fraction_non_alphanumeric': 0.042136915687116805, 'fraction_numerical': 0.01589409145936658, 'mean_word_length': 4.545257687310524, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Modal Approach to the Numerical Calculation of Primordial non-Gaussianities

We propose a new method to numerically calculate higher-order correlation functions of primordial fluctuations generated from any early-universe scenario. Our key-starting point is the realization that the tree-level In-In formalism is intrinsically separable. This enables us to use modal techniques to efficiently calculate and represent non-Gaussian shapes in a separable form well suited to data analysis. We prove the feasibility and the accuracy of our method by applying it to simple single-field inflationary models in which analytical results are available, and we perform non-trivial consistency checks like the verification of the single field consistency relation. We also point out that the i epsilon prescription is automatically taken into account in our method, preventing the need for ad-hoc tricks to implement it numerically.


Introduction
The deviation from perfect Gaussian statistics of the primordial cosmological fluctuations enables us to discriminate amongst the candidate physical mechanisms that produced the seed primordial fluctuations (for recent reviews, see for instance [1][2][3][4][5]). With the imminent arrival of Planck data together with those from large scale structures surveys, it is fundamental that theorists have firm predictions in wide classes of models to match against these increasingly precise data. However, most non-Gaussian shapes up to now have been determined under restrictive approximations, typically of the slow-varying type, to enable an analytical calculation. Interesting inflationary dynamics nonetheless exist which violate these assumptions and for which no theoretical prediction is yet available, for instance models with a strong background time-dependence or multifield models with non-trivial trajectories (see for example [6] for a recent illustration).
Procedures to numerically calculate the non-Gaussian properties of early-universe models do exist [7][8][9][10], but they rely on direct k-space configuration by k-space configuration calculations requiring tricks to implement the i prescription necessary to project onto the interacting vacuum, like the introduction of a 'damping' factor into the required time-integral to ensure early-time convergence.
In this paper, we propose an alternative method based on the separability of the treelevel In-In formalism and the modal decomposition techniques of Fergusson, Shellard and collaborators [11][12][13][14][15][16][17][18][19][20] (see also [21][22][23]). It enables one to numerically compute tree-level higher-order correlation functions of cosmological observables from any early-universe scenario, with three main advantages: • no ad-hoc trick is required to implement the i prescription, which is automatically taken into account in our method.
• the results of our procedure are smooth functions of their arguments, the k i 's, and not values of the relevant correlation functions on a grid of points in Fourier space.
• the resulting shapes are directly represented in a separable form, and hence can be efficiently constrained by data, for instance from the Planck satellite.
The outline of the paper is as follows: in section 2, we present the key-idea behind our method, namely the use of the intrinsic separability of the tree-level In-In formalism, and details about its numerical implementation. In section 3, we demonstrate its efficiency by calculating the tree-level bispectra generated by various operators in single-field test scenarios, like chaotic and Dirac-Born-Infeld inflation, and comparing them to known analytical results. We also perform non-trivial checks of our formalism, including the verification of the single field consistency relation between the scalar spectral index and the squeezed limit of the primordial bispectrum [24,25]. We state our conclusions in section 4.

Separability of the In-In formalism
The central starting point of our method relies on the realization that the In-In formalism used to calculate primordial cosmological correlation functions [26] is intrinsically separable, at least at tree-level.
The In-In (also known as Keldysh-Schwinger) formalism gives the following prescription for calculating the expectation value of any operator Q(t) in the interacting vacuum of the theory: where |0 denotes the vacuum of the free theory, T the time-ordering product, H I the interacting Hamiltonian, Q I (t) is in the interacting picture, and the i indicates that the time integration contour should be slightly rotated into the imaginary plane in order to project onto the interacting vacuum [24].
Concentrating, for simplicity of presentation, on the bispectrum of the curvature perturbation ζ generated in single-field models of inflation, one is led to evaluate, at tree-level: where τ denotes the conformal time such that dτ = dt/a(t), H (3) is the cubic interacting Hamiltonian, and all operators on the right-hand side are in the interacting picture (we omit the I for brevity). In practice, the dimensionless quantity we will compute and which is ultimately the most useful for data analysis is the so-called shape S(k 1 , k 2 , k 3 ), defined such that [12,27] where A s is a normalization factor related to the amplitude of scalar fluctuations. Illustrating this on the two cubic interactions where g(t) denotes a generic time-dependent coupling constant, we find respectively that where ζ k denotes the time derivative of the mode function with respect to τ , and k 2 · k 3 = 1 2 (k 2 1 − k 2 2 − k 2 3 ) using the triangle condition k 1 + k 2 + k 3 = 0. The key-point to notice is that, up to explicitly known k-dependent factors, like k 2 · k 3 in Eqs. (2.3)-(2.4), the primordial shape generated by any cubic operator is explicitly given in a separable form, i.e. as the time integral of a product of three functions of respectively k 1 , k 2 , k 3 . This crucial fact therefore enables one to apply efficient separable techniques to determine a modal decomposition of the non-Gaussian profile S generated by any cubic operator.

Modal decomposition
Following the modal techniques of Fergusson, Shellard et al, we wish to represent any primordial shape S(k 1 , k 2 , k 3 ) (or similar quantities related to the relevant shape by multiplication by known k-dependent factors and/or symmetrizations) as where the Q n 's are suitable basis functions that are separable and orthonormal, so that the α n 's can be computed efficiently. For this purpose, we should halt and define both the space over which we want to expand the shapes and the scalar product with respect to which orthonormality is defined. Physically, the shapes S are defined on a domain in Fourier space such that: (i) all k i 's ∈ [k min , k max ], where k min and k max stand for respectively the largest and smallest scales of interest (which depend on the particular probe we are interested in, for example CMB or large scale structures), (ii) the wavenumbers are such that the triangle condition k 1 + k 2 + k 3 = 0 holds, namely However, it is easy to realize that calculating the coefficients α n 's on such a tetrahedral domain would come at the cost of losing separability. Therefore, we consider instead the expressions (2.3)-(2.4) (and similarly for other vertices) as mathematically defining shapes over the cubic domain C ≡ [k min , k max ] 3 , irrespective of the triangle condition. Of course, for physical interpretations, the shapes we will thus construct should be restricted in the end to the tetrahedral domain, but the use of shapes defined on C enables us to use the intrinsic separability of the In-In formalism. Indeed, using the canonical scalar product on C: orthonormal and complete families of separable basis functions Q n with the correct symmetry properties can be constructed out of products of Legendre polynomials. For instance, a completely symmetric function of (k 1 , k 2 , k 3 ), like the building block of Eq. (2.4): can be represented on C as in Eq. (2.5), using the basis functions where P i denotes the i-th order Legendre polynomial, the index n represents the symmetric triplet {n 1 , n 2 , n 3 }, and where ∆ n 1 n 2 n 3 equals 1, √ 3, √ 6 if all the n i 's are equal, two of them are equal or all of them are different respectively. Using these orthonormal and separable basis functions, together with the separable representations of shapes as given by the In-In formalism, the modal decomposition coefficients α n 's in Eq. (2.5) can be computed as α n = CS Q n , where the triple integral over the cube C can be crucially decomposed as products of three one-dimensional integrals over k 1 , k 2 and k 3 respectively: and where the argument of the P n 's are (2k n − (k min + k max ))/(k max − k min ). Of course, a similar expansion holds for non-symmetric functions, like the building block of Eq. (2.3) which is symmetric only in (k 2 , k 3 ). The only difference in such cases is that the basis polynomials should only be symmetric in (k 2 , k 3 ). Following these methods, the shape generated by any cubic operator can be represented by a modal decomposition whose coefficients can easily be calculated numerically. Finally, let us note that, although we have illustrated our reasoning by considering the curvature perturbation ζ only, it should be clear that similar techniques based on separability hold in theories with several degrees of freedom, like in multiple field inflationary models.

Numerical considerations
It turns out that most of the difficulties we encounter to numerically evaluate the Fourier space integrals needed to determine the modal decomposition coefficients can be understood by using the crudest approximation for the mode functions, namely up to some slowly varying multiplicative factor function of both k and τ , and where c s is the slowly-varying speed of sound of fluctuations. Required integrals of the type kmax k min in (2.9) can be seem as sum of terms of the form kmax k min for n ≥ 0. An obvious difficulty in numerically evaluating these integrals is the rapid oscillations coming from the phase factor e −ikcsτ in the mode functions. However, this can be circumvented by numerically solving the linear mode function equation not for ζ k (τ ) itself, but rather forζ k (τ ) = ζ k (τ )e ikcsτ , which displays no rapid oscillations. This results in the integrand in Eq. (2.13) being put in the form of the slowly varying function k 2+nζ k (τ )ζ * k (τ ) multiplied by the explicitly known phase factor e ik(cs(τ )τ −cs(τ )τ ) , a type of integrand for which known efficient numerical routines can be used. In practice of course, the above slowly-varying function of k is only known as a result of evaluating the mode functions at sample k-points and constructing an interpolating function out of them. This results in one last difficulty: for n = 0, k 2+n ζ k (τ )ζ * k (τ ) contain terms behaving like 1/k which do not admit a rapidly convergent polynomial representation. We overcome this by simply considering the well behaved function k 3 ζ k (τ )ζ * k (τ ) (instead of k 2 ζ k (τ )ζ * k (τ ) in Eq. (2.12)), and dividing by the appropriate factor of k at the end of the calculations. On the other hand, note that no such trick is required for evaluating the integrals of the form kmax k min P i k 2 ζ k (τ )ζ * k (τ ) dk because of the additional k 2 factor in ζ k (τ ) in Eq. (2.11). Eventually, we observed that a slightly different expansion leads to improved accuracy, which consists in considering the shapes as function of the y i = log(k i )(i = 1, 2, 3), rather than the k i 's themselves, and therefore in decomposing them similarly as described above, but with the scalar product S(y 1 , y 2 , y 3 )S (y 1 , y 2 , y 3 ) dy 1 dy 2 dy 3 . (2.14) In the following, we denote results obtained by this method as resulting from the logmethod.

Projecting onto the interacting vacuum and the i prescription
When directly computing the value of a non-Gaussian shape in a specific momentum configuration, using e.g. Eq. (2.4), the i prescription is needed to project on to the interacting vacuum. However, the quantities we are calculating in our approach are different: these are the coefficients that enter into the shape's modal decomposition (2.5), like α n in Eq. (2.9). Here, the time-integrand comes as the result of integrating the mode functions in Fourier-space, which averages out the rapid oscillations of the mode functions deep inside the horizon, with one important consequence: evaluating Eq. (2.9) for > 0 and taking the limit → 0, following the i prescription, is equivalent to just evaluating it at = 0 straight away. Or, to put it a different way: the process of first integrating over Fourier space before performing the required time-integral is equivalent to the i prescription for us. Hence, while integrals like the one in Eq. (2.4) do not converge, as their lower bound τ ini goes to −∞, without an appropriate rotation of the time integration contour into the imaginary plane, this is not the case in our method. Actually, by using the asymptotic form of the mode functions (2.11) deep inside the horizon -which only results from our choice of the Bunch-Davies vacuum -and the behavior of the scale factor a(τ ) ∝ 1/τ , one can easily determine the early-time behavior of the required time-integrands, like the one in Eq. (2.9). One thus finds that the coefficients of the modal decomposition converge, as τ ini goes to −∞, like ∞ dx x 2 e ix for the operators in dt a(t) ζ 3 and dt a 3 (t) ζζ 2 (whatever their explicit gradient structures), and like ∞ dx x e ix for the operators in dt a 3 (t)ζ 3 and dt a(t)ζζ 2 . In practice, these behaviors, which we have verified in our numerical calculations, imply that requisite time-integrations can be performed only starting from a few e-folds before the largest scale of interest crosses the horizon.

Results
In this section, we demonstrate that the method described above works as intended. We first describe the details of our set-up, in particular the inflationary backgrounds and the cubic operators we have considered, before showing explicit results in various situations.

Set-up
Test scenarios.-We consider two kinds of inflationary scenarios: the one of a canonical scalar field with an m 2 φ 2 potential, which we refer to as the chaotic scenario, and the one corresponding to a Dirac-Born-Infeld (DBI) action in an AdS geometry [28]. In more detail, the two corresponding actions read (we omit the Einstein-Hilbert action which is common to the two set-up, and use units in which M p = 1) Our aim is not to consider inflationary scenarios that match precisely observational data, for example in terms of the amplitude and spectral index of the scalar fluctuations they generate or their level of non-Gaussianities, but rather to consider backgrounds in which all the slowly-varying approximations that are usually made in analytical treatments are satisfied, so that our numerical procedures can be tested against reliable theoretical predictions.
Unless otherwise stated, we take m = 10 −5 in both scenarios and λ = 10 16 . We also use initial conditions such that background attractor solutions are reached rapidly, namelẏ As for the linear fluctuations about these inflationary backgrounds, the curvature perturbation obeys the simple second-order equation where it is understood that c 2 s ≡ 1 in the chaotic case, and we impose initial conditions such that the Bunch-Davies behavior holds inside the horizon: Eventually, we consider a range of scales such that k min c si a i H i = 10 3 -all studied scales are therefore well inside the horizon initially 1 -and with k max /k min equals to 10 2 (similar results are obtained with k max /k min = 10 3 ).
Cubic operators.-Using the background and the linear fluctuations just described, we would like to test our method for numerically calculating the non-Gaussian shapes generated by various cubic operators. As we have explained in section (2), we factor out any explicitly known k-dependence coming from gradient terms in cubic operators. We therefore have to apply our method only to four kinds of cubic combinations of ζ, namely ζ 3 , ζ 2ζ , ζζ 2 andζ 3 , from which one can then reconstruct the shapes of any specific operator. In the context of the most general second-order scalar-tensor theory, also known as Horndeski theory [29][30][31], it has been shown in Ref. [32] that only 5 independent cubic interactions can possibly be generated, i.e. that the third-order action can effectively be put in the form where the g i 's are time-dependent coupling constants, and with other interactions, like for example dt d 3 x aζ(∂ζ) 2 , being redundant, in the sense that they can be expressed in terms of these 5 by using the linear equation of motion (3.3). This fact provides the opportunity of interesting consistency checks, and we have therefore additionally considered this latter operator.
Ordering of the basis functions.-We have described in subsection (2.2) our choice for the basis functions Q n we are using to expand any shape in a separable form. Of course, in any practical application, we should truncate the infinite expansion (2.5) and use a specific ordering for the Q n 's. For that purpose, we use the 'slicing' ordering described in Ref. [13] and decompose functions with polynomials (of k i or log(k i ) depending on the chosen method) up to a specific total order N max . For symmetric functions for example, the index n thus relates to the completely symmetric triplets {n 1 , n 2 , n 3 } (with a specific choice of sub-ordering): where the transitions between total polynomial order is underlined.

Numerical results
Bare operators.-To show the effectiveness of our method, we applied it to the six operators listed above, with constant g i 's in Eq. (3.5), for the two different test scenarios and using N max = 6 in the modal decomposition 2 . We then evaluated the correlationsas measured by F (S, S )/ F (S, S)F (S , S ) with the scalar product F (2.7) restricted to the tetrahedral domain -between the corresponding theoretical and numerically calculated shapes. Our results are listed in table (1). The various cosines are impressively close to 1, demonstrating the ability of our method to reconstruct non-Gaussian shapes. operator correlation (m 2 φ 2 ) correlation (DBI) ζ 3 0.9992 0.9994 ζ(∂ζ) 2 0.99997 0.99995 0.999990 0.99998 Table 1. Correlations between theoretical and numerically calculated shapes generated by various (bare) operators on our two test background inflationary scenarios.
Additionally, its rapid convergence is displayed in table (2), where we indicate the correlations between theoretical and numerically calculated shapes for increasing N max . We give the results for two representative operators of the equilateral and local type respectively, namelyζ 3 and ζ(∂ζ) 2 , and for the test chaotic inflationary background (similar results are obtained in the DBI case). Eventually, let us note that it takes about 1.5 minutes on a laptop to calculate all possible bispectrum shapes (remember that they can be all deduced from the 4 cubic combinations of ζ andζ) with N max = 2, and 4 minutes with N max = 6. While this is sufficient for our purpose, we expect that a significant speed-up can be obtained after optimization of our code. Two non-trivial checks.-In single-field inflation, there exists a well known consistency relation between the amplitude of the primordial bispectrum in the squeezed limit and the scalar spectral index [24,25], namely that where n s (k) is defined as To demonstrate the accuracy of our method, we calculated the bispectrum generated in various chaotic inflationary scenarios, which, at leading-order in a slow-varying approximation, simply corresponds to setting g 2 (t) = ( − η) and g 3 (t) = ( + η) in Eq. (3.5), where η ≡˙ /(H ). We then evaluated the reconstructed bispectra on squeezed triangles such that k 3 = 0.01k 1 = 0.01k 2 and we numerically calculated the corresponding n s (k 1 ). We found that in all cases the single field consistency relation was verified to better than the % level. For instance, for φ i = 16, we find that n s (k max )−1 = −0.0379, to be compared with −0.0377 deduced from the numerically calculated bispectrum, so a 0.5% discrepancy.
As we mentioned above, another non-trivial check of our numerical method involves the redundancy of some cubic operators [32]. At leading order in the slow-varying approximation, the DBI third-order action reads [33] (3.10) Surprisingly enough at first sight, the two cubic operators present in this action generate non-Gaussian shapes that are strongly correlated with the local ansatz [34], whereas the DBI scenario is known to generate a bispectrum of the equilateral type [35]. A way of understanding this is to note that, upon using the linear mode function equation (3.3), an equivalent DBI third-order action reads, at leading-order in the slow-varying approximation 3 [32]: This latter form of the third-order action makes it manifest that the total DBI bispectrum is of equilateral type, as this is the case for the individual bispectra generated by the operators inζ 3 andζ(∂ζ) 2 . The fact that the local-type behaviors of the two operators ζζ 2 and ζ(∂ζ) 2 should conspire to cancel in (3.10), resulting in a total equilateral-type shape, provides a nontrivial test of our method. We have verified it, obtaining a correlation of 0.993 between the bispectra calculated using Eq. (3.10) and Eq. (3.11) in our test DBI scenario. This is summarized graphically in Fig. (1), where shapes are represented through several density contours.
Other examples.-For illustrative purposes, we display below examples of theoretical bispectrum shapes (left) and their numerically calculated counterparts (right) in our test DBI Figure 1. A graphical representation of a non-trivial check of our method. From left to right: the shapes generated in our test DBI scenario by, respectively, the ζζ 2 and ζ(∂ζ) 2 operators in Eq. (3.10), the resulting total shape, and the total shape calculated from Eq. (3.11). Each shape are represented through several density contours. The third one shows that the cancellations of the local-type behaviors of the first two shapes is not perfect in the squeezed limits (in blue). The overlap with the fourth shape is nonetheless remarkable, with a correlation of 0.993 between the two. scenario. We used the log-method (c.f. subsection (2.3)) with N max = 8, together with the coupling constants of Eqs. (3.10)-(3.11). The agreement is impressive.

Theoretical bispectra
Numerically calculated bispectra Note that because the slow-varying approximation is very accurate on our test DBI background, the numerically calculated shapes displayed above are almost exactly scale invariant. For comparison, we show in Fig. (2) the shapes generated by the two operatorṡ ζ 3 andζ(∂ζ) 2 in the third-order action Eq. (3.11) under the same conditions as previously, but with λ = 10 12 in the action (3.2). This modification results in a much more rapidly varying background, and hence to strongly scale-dependent equilateral non-Gaussianities, whose amplitude in the equilateral limit grows by a factor of about 30 between the largest and the smallest scales. This exemplifies the effectiveness of our approach for studying primordial non-Gaussianities generated from non-trivial inflationary dynamics. Figure 2. Shapes generated by the two operatorsζ 3 (left) andζ(∂ζ) 2 (right) in the third-order action Eq. (3.11), on the DBI background with λ = 10 12 (see the main text). The rapidly varying background leads to strongly scale-dependent equilateral-type bispectra.

Conclusions
In this paper we have proposed a new method to numerically calculate higher-order correlation functions of primordial fluctuations generated from early-universe scenarios. It is based on the realization that the tree-level In-In formalism is intrinsically separable, which enables one to use modal techniques to efficiently represent any non-Gaussian profile. We have proved the feasibility and the accuracy of our method by applying it to simple single-field inflationary scenarios in which analytical results are available, finding impressive agreements and rapid convergence with respect to the number of modes used in the decomposition. We also performed non-trivial consistency checks, including the verification of the single field consistency relation, and the cancellations of two local-type behaviors to obtain an equilateral-type shape, as relevant in the DBI inflationary model. Additionally, we found our code to be quick, although a significant speed-up can be expected after proper optimization.
To us, the main advantages of our method are that: (i) the i prescription is automatically taken into account, preventing the need for ad-hoc tricks to implement it. (ii) We are able to calculate correlation functions at once over the whole domain of interest, and not only for specific k-space configuration. (iii) Non-Gaussian shapes are directly expressed in a separable form, and hence are ready to be efficiently constrained by data.
While we have mainly applied our method to simple early-universe models, thus providing a proof of its concept, its usefulness ultimately lies in the fact that it provides an efficient way of determining the non-Gaussian properties of cosmological fluctuations in non-trivial models in which analytical techniques are known to be insufficient. The result of such investigations will be presented elsewhere [36,37].