On the initial condition of inflationary fluctuations

It is usually assumed that the inflationary fluctuations start from the Bunch-Davies (BD) vacuum and the $i\varepsilon$ prescription is used when interactions are calculated. We show that those assumptions can be verified explicitly by calculating the loop corrections to the inflationary two-point and three-point correlation functions. Those loop corrections can be resumed to exponential factors, which suppress non-BD coefficients and behave as the $i\varepsilon$ factor for the case of the BD initial condition. A new technique of loop chain diagram resummation is developed for this purpose. For the non-BD initial conditions which is setup at finite time and has not fully decayed, explicit correction to the two-point and three-point correlation functions are calculated. Especially, non-Gaussianity in the folded limit is regularized due to the interactions.


Introduction and Summary
Inflation is the leading paradigm of the early universe cosmology. The fluctuations generated during inflation provide seeds for the cosmic microwave background (CMB) and the large scale structure (LSS) formation [1] . The standard calculation of those primordial fluctuations follows from the quantum theoretical in-in formalism following two assumptions, namely the standard vacuum initial condition (known as the Bunch-Davies vacuum [2], or BD vacuum for short) and the iε prescription. Those assumptions are inherited from the flat space quantum field theory, but have to be reconsidered in cosmology.
• The BD vacuum initial condition. This is the simplest choice of initial state in the simplest models of inflation, because inflation is an attractor solution. However, it has been debated for long because of the following issues: -Beyond the attractor stage of inflation. Scale dependent features during inflation can override the BD vacuum initial conditions . For example, inflation may be just enough and the start of observable inflation may be close to the absolute start of inflation [3]. As another example, there may be features on the inflationary potential such that the inflationary fluctuations after the features are in an excited state before horizon crossing [4].
-Beyond the simplest theory of fluctuations. In cosmology it is convenient to follow the time evolution of a comoving perturbation mode. The comoving mode originates from scales much smaller than the inflationary Hubble scale, and its physical wavelength expands with the cosmological expansion. During the expansion of its physical wavelength, the dynamics of the mode may be governed by different effective field theories, or no effective field theory at all when its physical wavelength is shorter than the Planck scale (the trans-Planckian problem [5]).
-Beyond the leading order calculation of gravitational fluctuations. Practically, the BD vacuum is selected as the lowest energy state. However, when gravitational fluctuations are concerned, energy is a gauge dependent quantity. Different gauge can have different definition of time, and thus different definition of energy. This is similar to the case that in the Minkowski vacuum, accelerating observer sees Unruh radiation [6], which appears to be no longer the lowest energy state. In the literature in each gauge people choose the lowest energy state as the physical "vacuum" state. This cannot be right. Only the vacuum state of one gauge should be physical and the vacuum in other gauges should be the gauge transformation of the same physical vacuum.
• The iε prescription. This prescription is not relevant in the tree level power spectrum calculation, but become important for the non-trivial in-in calculation for higher point correlation functions or loop diagrams. In flat space in-out formalism, the iε prescription is proposed to project the physical interacting vacuum onto the vacuum of the free theory, because only the vacuum of the free theory is operationally defined by the free quantum fields (or the interacting picture fields) and can be practically used in the perturbative calculation. One can relate the free vacuum |0 and the interacting vacuum |Ω by where T is the duration of the interaction, H is the full Hamiltonian, E 0 is the energy of the ground state defined by E 0 ≡ Ω|H|Ω , and E n ≡ n|H|n for non-perturbative states |n with higher energies. One can then send T to ∞ by T → ∞(1 − iε). Then all but the first term in the RHS of (1.1) vanishes, and we obtain a relation between |0 and |Ω . The following assumptions are involved in this prescription: -One can adiabatically turn off the interactions. This assumption works fine in flat space calculation of the S-matrix because we are preparing the initial states in the far past with large spatial separation. Following cluster decomposition [7], or any explicit law of forces, the states can indeed be considered to be non-interacting.
Actually, under some mild assumptions, the validness of relating interacting vacuum to the free vacuum in this way can be rigorously proved in quantum field theory, known as Gell-Mann and Low theorem [8]. However, in cosmology, we are interested in considering the time evolution of the initial vacuum state. The state is initially of sub-Hubble size and all (virtual) particles stay close to each other. Thus we are no longer sure about the validity of turning off interactions in the calculation of cosmological perturbations.
-There exists enough time duration T for the iε prescription. This assumption is again tricky in cosmology, because this statement is again coordinate dependent. For inflation, one can use conformal time or proper time. When the conformal time is used, one indeed have nearly infinite (though still not really infinite because inflation cannot be eternal to the past) amount of conformal time in the past. However, when using proper time, the amount of time duration gets shortened exponentially. One can indeed argue that before horizon crossing, the conformal time is more relevant. But explicit calculation is needed to verify the argument. Even we use the conformal time, a mathematically infinitesimal i does not work because of the finiteness of conformal time, even if the duration is exponentially long.
In this work, we aim to provide a systematic method towards resolving the above puzzles. This is an extension of our previous work [9]. We show that interaction is the key to the vacuum and the i problems.
Interaction exists in the early universe. The theory of gravity is nonlinear. The gravitational nonlinearity provides a lower bound on the interaction of perturbations during inflation. In terms of the non-Gaussianity estimator f NL , the minimal gravitational nonlinearity corresponds to f NL ∼ O(0.01). Large non-Gaussianities are predicted in some inflation models and the current observational bound is f NL of order 10 or 100, depending on the shapes of non-Gaussianity.
For this purpose, in our previous work, we calculate the one loop correction of the two point function with non-BD initial conditions. We have shown that, with the help of interactions, the non-BD initial conditions dissipates exponentially fast towards large scales. The one loop correction of the non-BD coefficients can be classified into two types, namely the correction to the amplitude and phase of the non-BD coefficient. The correction to the amplitude of non-BD coefficients corresponds to the contributions close to the folded limit of the interaction vertex. This amplitude correction is negative and can be resumed onto the exponent by dynamical RG method [10][11][12]. As a result, for sub-horizon fluctuations, we have where c k is the absolute value of the tree level non-BD coefficient, c eff k is that with dynamical-RG-resumed one loop corrections, and τ 0 is the initial time where the non-BD initial condition is setup. For f NL ∼ O(1), the characteristic scale on the exponent is between kτ 0 ∼ 4, indicating that non-BD initial conditions which are setup at sub-Horizon scales as deep as 4 e-folds start to decay exponentially. For larger non-Gaussianities, the decay of non-BD initial conditions become significantly faster. As a result, smaller non-Gaussianities, which seem not great for the purpose of probing interactions during inflation, have the advantage of better preserving the initial state of inflation.
In this work, we solidify the previous calculation by an explicit loop calculation, fixing the previously undetermined order one coefficient. For (∂ t ζ) 3 interaction, the result is The dynamical RG resummation method which has been used in our previous work is also checked explicitly using a direct resummation of one particle reducible multi-loop diagrams. We show that the two results agree up to a two-loop contribution, which is under control when proper scale of renormalization is chosen. There are model dependent and model independent components in (1.3). The numerical factor is of course model dependent. The 5th power in τ is also model dependent. If the interaction were marginal (i.e. dimension 4 after canonically normalize ζ), then one expects linear dependence in τ , because the total amount of interaction should be proportional to the length of interaction time. Here, the operator under our consideration has dimension 6. Thus for each interaction vertex there arises two additional powers of τ due to UV sensitivity. As a result the exponent scales as τ 5 . For inflation with standard kinetic term and Einstein gravity, the interactions have dimension 5 and we should expect the exponent scaling as τ 3 . The dependence on P ζ , f NL and the exponential structure of the decay, on the other hand, should be model independent. Also, the interaction scales linearly in τ − τ 0 when τ − τ 0 is small. This is model independent from the physical interpretation of a decay rate.
Technically, it is interesting to note that, in the sub-horizon limit, the reducible multiloop diagrams (as a chain of one loop diagrams) dominate over the irreducible ones. The reason is as follows. We hope to pick up the highest power of |kτ 0 | in the calculation. The highest power comes from the diagrams where the largest number of vertices can freely take values from τ ∼ τ 0 to |kτ |∼ 1, which is a large range. In the reducible diagrams, the vertices group into freely moving pairs, each pair represent a loop and the relative time difference is constrained by the uncertainty principle. However, for diagrams which contain irreducible multi-loop parts, more vertices are constrained by the uncertainty principle and thus do not show up at leading power of |kτ 0 |. This further assures the validity of the dynamical RG method.
We then study the one loop correction of the three point correlation function. In the case of the three point function, the three external legs can carry different momenta and thus the dynamical RG method becomes no longer accurate. We can nevertheless still calculate the multi-loop reducible diagrams and sum them up explicitly. The result corresponds to adding a decaying factor to the propagator: With the help of the resumed propagator, the folded limit of non-Gaussianity no longer diverges. The folded contribution of non-Gaussianity vanishes if taking τ 0 → −∞. Once a finite initial time τ 0 is given, explicit loop-corrected shapes of non-Gaussianities can be obtained. For example, if the non-BD modes are set up at relatively early time, the non-Gaussianity may show some nontrivial shape like Figure 1. One of the underlying reason is that the large k modes decay faster, while small k modes decay relatively slowly and thus leave more prominent non-BD initial information on the observations. Also, it is known that once given an initial time τ 0 , one can no longer use the iε terms to suppress the boundary terms in the UV. As a result, the tree level result depends on τ 0 strongly and oscillations are present if the cutoff is sharp. Such dependence presents also for the BD initial condition. We show here that those τ 0 dependent terms also decay with a similar exponent. In other words, the interactions practically serve as the iε, and indeed pick up the physical initial state.
The rest of the paper is organized as follows: In Section 2, we write down a simple interaction model and review the basic formalism. In Section 3, we calculate the two point correlation function. After recapitulating the one-loop folded limit cut-off result as given in our previous work, we come up with the new technique: loop chain diagram resummation. The dynamical RG method is also used to double check the result. In Section 4, we calculate the loop corrected three point function. The decay of both the non-BD terms and the non-iε suppressed terms are manifest.

Our Model
We start from general single field inflation with L = P (φ, X) [13,14]. The second and third order action up to the first order in slow parameter can be derived as where the dot "˙" denotes the derivative with respect to time t and the prime " " denotes the derivative with respect to conformal time τ . Also we set the reduced Planck mass M p = 1.
The relevant quantities of this model are From non-interacting S 2 , we can quantize the field ζ: with superscript "I" for interacting picture. The mode function is given by As our motivation is to see the effects of interactions, we can use the sub-horizon limit approximation |c s kτ | 1. The reason is that after horizon crossing, the modes are nearly frozen and can not evolve anymore. So, interactions can play no role in the super-horizon case.
In such limit, the mode function and its derivative have the following approximate behaviors u k ∝ kc s τ e ±ikcsτ ,u k ∝ 1 a k 2 c 2 s τ e ±ikcsτ . (2.6) For ζ field, the leading order time dependence is similar and note that the real space derivative corresponds to momentum multiplication in momentum space ∂ζ ↔ kζ k . So, we have the following relation 1 In the 3rd order action S 3 for interaction, the ratios of different terms are 2nd term 3rd term (2.9) We can easily see that in the sub-horizon limit, usually the first term is much larger than other two terms. So, one can just consider the first term and discard other two. This is because the highest dimensional operator is the most sensitive to the UV physics. Furthermore, for simplicity, we set sound speed c s to be 1. Based on these arguments, we can consider a simple model of inflation described by The Hamiltonian for interaction is where the scale factor a(τ ) ≈ − 1 Hτ for quasi de-Sitter space during inflaton. In this simplified model, the mode function is given by where the coefficients C + , C − are subject to the following constraint required by the consitency of quantization In the usual case, the requirement of the vacuum state as a minimal energy state or the matching of de-Sitter space-time in the sub-horizon limit with Minkowski space-time will give rise to another condition C − = 0. This is the so called Bunch-Davies vacuum [2]. But, here we consider small C − , corresponding to non-Bunch-Davies case. To the first order, we have C + (k) ≈ 1, C − (k) ≈ c k e iθ k . So the mode function and its derivative are given by (2. 16) In order to make the story simple and clear in some sense, in the following calculations, we assume that the mode functions do not depend on the directions of momentum. Namely, we require u k = u k , θ k = θ k . The calculations and conclusions are expected to be more general independent of these assumptions except possible complications.
3 Two-point function

General consideration: tree level and one-loop level
The interaction Hamiltonian in interaction picture is (3.1) The two-point correlation function can be calculated by using the in-in formalism (see Appendix A): The zeroth order of two-point correlation function is given by (note our notation k = k 1 ): where we consider the sub-horizon limit −kτ 1 and only keep terms up to the first order in c k .
The first order correction of two-point correlation function vanishes due to odd number of operators or imbalance of creation and annihilation creators.
Next, we consider the second order loop corrections. There are two types of corrections: non-BD mode in the loop and non-BD mode in the external line. When the non-BD modes are in the external line, the physical meaning is very obvious if we cut the loop. This process can be thought as the decay of non-BD mode in the external line into two BD modes in the loop. Furthermore, in order to match with tree level result as will be elaborated later 1 , we need to have something like c k which also implies a non-BD mode in the external leg. Thus, for simplicity, we can just consider this case by setting c p , c q = 0.
We are interested in the sub-horizon limit which means that |kτ | 1. Usually, this doesn't imply |(p + q − k)τ | 1 in the folded limit. But, for simplicity, let's first consider the unfolded case. In such a case, when evaluating the above equations, we only keep those terms which have highest power in τ and zeroth and first order in c k . Then we can use the following integration formula τ n e iQτ dτ ≈ τ n iQ e iQτ + · · ·. The second order symmetric part is: The second order asymmetric part is: The physical meaning of matching the tree level result is that, the contributions coming from cp and cq correspond to processes where two long modes fuse into a short mode. When the short mode is far away from vacuum and the long mode is nearer to the vacuum (considering more time of decay), this is unlikely to happen. However, there is an important exception: Near thermal equilibrium, the detailed balance makes sure that the decay of the short mode is indeed balanced by the fusion of the long mode. Our calculation thus does not apply for such near equilibrium cases. An approach of Boltzmann equation would help and we hope to explore this possibility in the future.
where we use the approximation u k (τ ) It should be noted that we have τ 5 terms now. This is because the exponential parts of mode functions cancel and the power of τ increases after integration. So, the final result is Because we are considering the sub-horizon limit, τ 5 terms dominate. Note that in Eq. (3.7), the oscillation is sine function form while the tree level result in Eq. (3.3) is cosine function form. This means that the unfolded part can not modify the amplitude of effective non-BD coefficient at leading order.
With this in mind, we need to reexamine the asymmetric part: where the ellipsis denotes the relatively irrelevant terms. Now that the tree level result has different function form with unfolded loop correction, we may expect the dominant loop correction comes from the folded limit. So, we can Taylor expand the integrand and perform time integration, yielding (3.10) As we expect, the tree level function form emerges which would be the dominant loop contributions. After that, loop momentum integration can be performed as follows by choosing a momentum cut-off near the folded limit (p + q − k) ≤ Λ, When evaluating the above momentum integral, we use the elliptical coordinate system (see Appendix B). While for the unfolded part, we can use Eq. (3.7), so p+q−k>Λ . Naively, the coefficient is infinity due to non-bounded momentum integration. But physically, if we take the renormalization counter term into considerations, the coefficient should be finite. Nevertheless, its dependence on non-BD coefficients are different from the folded limit one and tree level result.
Note that in order for the expansion to be valid, we require spite of an order one discrepancy. Collecting all the facts, we get final two-point correlation function under one-loop correlation, where the ellipsis includes the loop corrections of BD modes and higher order corrections to non-BD modes. It is very interesting to see that the loop corrections to the non-BD coefficients are negative, implying the decay of non-BD modes.

Rigorous treatment of momentum integral
From previous calculations, in the sub-horizon limit, by power counting, the τ 5 terms are the dominant one. From Eq. (3.8), after loop correction, In our sub-horizon limit approximations, f S + f A actually is given by Eq. (3.9). We can regard the right hand side as non-BD contribution with effective parameters running with time, so 14) The effective one can be written as where δc k = c eff k −c k , δθ k = θ eff k −θ k are very small when τ and τ 0 are very close. By matching the form, we can get The integral I can be simplified by using elliptical coordinate system (see Appendix B) Let's focus on the integral I − first. We define z = µ − 1 and find that Then, we need to integrate over µ or z. From the expression for S(z), we can write it as the following general form The polynomial part of S(z) implies the power divergence if we integrate over z. They can be discarded if we choose to believe that these terms can be canceled by the the local counter terms. So, we only need to care the fractional part S f . In the UV region and sub-horizon limit, intuitively, the most relevant contributions come from those terms with lowest power in 1/z and highest power in τ , i.e.
for the current problem. But, this navie choice is a little problematic due to the divergence near z ∼ 0. Fortunately, we can simplify this problem by introducing a new basis functions T n (z) (see Appendix C) We can easily solve the above equation to obtain A n . In particular, A 1 is given by Thus we have A 1 = b 1 as emphasized in Appendix C. So, As mentioned in Appendix C, T 1 is not integrable and thus, we need to choose a momentum cut off on z, In the sub-horizon limit, we only need to keep the highest power term in τ , i.e. τ 5 terms. Then, I − can be approximated as It's interesting to note that Similarly, we can get Compared with the folded limit cut-off result (3.12), they only differ by a decay factor, i.e.
Another important observation is that τ 5 terms appear in unfolded case, folded limit cut-off result and present rigorous treatment with universal power. So, our strategy can be like this: analyze the unfolded case first, extract the highest power terms and then transform back to the folded limit. Very amazingly, the transformations can be simplified by just changing 1/z to T 1 (z) with exactly the same coefficients (remember b 1 = A 1 ).
The one loop correction is small as long as the initial time and final time are close enough. Once the time difference becomes large, the one loop perturbation is not valid any more. We need to try to cure the secular growth with time, either through the dynamical renormalization group method or by turning to higher order loop analysis.

Dynamical Renormalization Group method
The one loop corrections of effective parameters are Once we come to realize that c k should be c eff k and running, the c k above should be replaced with c eff k (τ ) . This yields The idea essentially is the dynamical renormalization group method [10][11][12]. The effects of early time modes to later time modes through loop corrections can be viewed as the modifications of effective parameters, yielding the running effective parameters with time. For DRG method, the physical picture is very clear and enlightened. In the next section, we are going to provide another way which is more rigorous in mathematics.

Multi-loop analysis and loop chain diagram resummation
In the following, we are going to consider the higher loop corrections to the non-BD coefficients. In principle, there are infinite ways to draw the corresponding Feynman diagrams, nested or non-nested. As we see before, all the modes which run in the loop are BD modes. Non-BD modes in the loop will not affect the effective non-BD coefficients c k and thus are not considered. Furthermore, we only consider the non-nested loop-chain diagram which consists of loops connected in series. Due to the time sequence of interacting vertices, there are still lots of diagrams which have the loop-chain topology but differ in the time ordering. Thanks to the sub-horizon limit, we only need to keep the highest power term in τ . In such a limit, we are going to show that only those V-shaped diagrams dominate.
Let us analyze the basic component of the diagram first. For each loop, there are three possibilities of time sequence as shown in Figure 2.
In Figure 2(a), suppose τ U > τ 1 > τ 2 > τ L (τ U , τ L are upper and lower limit of integration and τ 1 , τ 2 are time of two interacting vertices of the loop) and note that loop modes are BD-modes, we have the following integral where the τ 5 terms are contributed by e i(p+q−k)∆τ , while O(τ 4 ) are contributed by other terms. In the sub-horizon limit, only τ 5 terms are most relevant, so in the original expression, we only need to keep e i(p+q−k)∆τ like terms. In particular, in such case, all the related modes are BD modes. In Figure 2(b), suppose τ U > τ 1 > τ 2 > τ L , similarly we get Similarly τ 5 terms are contributed by e i(p+q±k)∆τ . Now, we have two such terms. The first one corresponds to the folded limit case, while the second one vanishes when k(τ U − τ L ) O(1) according to the previous one-loop analysis. Thus, we only need to keep e i(p+q−k)∆τ terms and now, the late time τ 2 external mode is non-BD mode. For Figure 2(c), the analysis is nearly identical to that in Figure 2(b). Besides the behavior of building components mentioned above, there are several other interesting properties for the whole diagram: • The time sequence of these interacting vertices should be all time-ordered or anti-timeordered. There is no mixing. In another word, in the in-in formalism, the relevant contribution comes from ζ 2 H...H and H...Hζ 2 which are related by some complex conjugation.
The reason is that when we go along the loop chain, we start from time τ and finally go back to time τ . So, there must be some turning and extremal loops. These loops should be non-BD loops if we require the τ 5 terms. While we only keep terms up to first order in non-BD coefficients c k , there can only be one non-BD loop and the only possible configurations are V-shaped with all time-ordered or anti-time-ordered interaction time sequence.
• The two time of vertex at each loop should be consecutive. The reason is not to disturb the time integration and get as high power as possible.
All these conditions are verified by explicit calculations. Remember that we only keep the lowest order term in c k , which means that we can only have one non-BD mode in whole diagram. In order to get as high order terms as possible, the diagram can only be composed of lots of BD loops (as shown in Figure 2(a)) and one non-BD loop (as shown in Figure 2(b)). The final dominant diagram is V-shaped loop chain diagram (as shown in Figure 2(d)) where the tip of V is the non-BD loop (see Figure 2(b)). The non-BD mode has the earliest time and connects the non-BD loop with other BD loop.
So, the final contribution is of the form (The time ordering is Here, we make an assumption or approximation: Mathematically, it will give rise to very easy and interesting result. The other part starts from two-loop. We drop this two-loop contribution for the reason that shall be explained at the end of this subsection. As we calculated before (I 0 = 1), (3.38) By mathematical induction, we can show that under the previous approximation (3.39) Proof : Suppose it holds for I L−1 , then where the momentum integral can be obtained from the following: where we substitute z = µ − 1 after integrating ν and + · · · represents all possible terms which can not be be written as 1/(p + q − k) or 1/(µ − 1) or 1/z. Using the basis function in Appendix C and note that b 1 = A 1 , we have  there are three mode functions at each interacting vertex which are symmetric, implying a factor (3 × 3 × 2) L . Finally, the coupling factor of interaction and the numeric factor of mode functions should also be included. After taking all of these into considerations, we are led to a factor (3.43) So, the final result for the effective non-BD coefficient is Remarks on the approximation: The real part of one loop result (3.26) is exact. The calculations above show that higher loop corrections can be thought as the power of one loop result in some sense. But the even power of the imaginary part of one loop is real, which means that the imaginary part of one loop result is very important in higher loops. But anyway, we can regard them as the phase factor corrections e iγ with γ ∼ O(τ 5 0 , τ 5 ). The rigorous treatment of this phase factor is beyond our capability due to the log term in Eq. (3.24) which may require involved and subtle regularization and renormalization procedure [15,16]. But we can justify the log term from a physical perspective.
In one loop calculation Eq. (3.25), the cut-off divergence related part is Ci(kΛ(τ − τ 0 )) − log(kΛ(τ − τ 0 )). When τ − τ 0 → 0, it vanishes. Otherwise, the cosine integral function contribute nothing and we can just consider the log part. The divergence is expected to be canceled by the counter term which has the form .... log(Λ phy /µ) where µ is the physical renormalization scale and the coefficient is exactly as that of our loop result if the divergence is indeed canceled. Note that Λ phy = Λk/a(τ mid ) with τ 0 < τ mid < τ . So, the final result should be of the form − log(Λ phy a(τ mid )(τ − τ 0 )) + log(Λ phy /µ) = log( H µ τ −τ 0 −τ mid ). In order for the one loop level perturbation result to be valid, k 5 (τ 5 − τ 5 0 ) can not be too large and thus τ /τ 0 ∼ O(1) in the sub-horizon limit |kτ |, |kτ 0 | 1. After realizing this fact, the log term can only contribute finitely. So, the imaginary part of Eq. (3.25) is Im . After applying the arguments to each loop of multi-loop calculation, we can conclude that the final contribution of one loop imaginary part is a phase factor correction e iγ with γ ∼ O(τ 5 0 , τ 5 ). The exact form or coefficient is complicated, for simplicity, we ignore the phase factor correction and only consider the amplitude suppression.

Tree-level result
Following the standard in-in formalism, we can easily obtain the tree level three-point function (4.1)

Loop correction
It is well known that when the non-BD initial condition is assumed, the folded limit non-Gaussianity will blow up. From the two point function calculation, we know that the effective non-BD coefficient will decay with time due to the loop correction. We expect that when we include loop corrections to 3pt function, the divergence will be cured.
For the three point functions, there are three external momenta. If we only consider the loop chain diagrams which modify three legs separately, namely those diagrams where three legs only meet at the original tree level interacting vertex with time τ V , things will become very easy. Other diagrams, either nested or connecting different legs, are thought to only result in small corrections because their contributions are lower orders in τ , which can  For loop chain diagrams, we need to consider loop corrections on each leg. This means that even the BD three point functions will be corrected by loops. We first consider this case because it is easier to deal with due to the symmetry of three legs.

Loop correction to BD three point function
Explicit case studies show that the dominant diagrams have the following properties: where we use the tricks presented in the Appendix A and τ m = τ V is the 3pt-tree level interacting vertex.
• All loop vertices time are later than three point function interacting vertex τ V < τ L . This implies n < = 0.
• For each leg, the loop chain has similar properties as those stated before including consecutive loop time.
So, the loop corrected three point function can be written as For diagram with L 1 , L 2 , L 3 loops at each leg respectively, the result contains the following term L . Next, we need to include coupling constants and do combinatoric counting. For each leg, the result is So, the loop contribution factor with specific number of loops is (4.6) And finally, we need to sum over all possible loops, yielding This is just the philosophy presented in Figure 3, i.e. replacing "1" with some exponential suppression factor and then you get the loop corrected results.

Loop correction to non-BD three point function
If there is one non-BD mode in the diagram, things are very similar to BD one but a little more complicated due to different behaviour of non-BD mode.
For illustration, we consider the one-loop case first. This loop is in the non-BD leg. If we denote the loop time as τ 1 , τ 2 (τ 1 > τ 2 ), there are two diagrams with τ 1,2 > τ V or τ 1,2 < τ V . Diagrams with τ 1 > τ V > τ 2 means that the loop time are not consecutive and can not contribute highest power terms from previous analysis.
For three point function, there are two types of diagrams: • Type 1 ( Figure 5(b)): If τ V < τ 2 < τ 1 < τ , the non-BD mode is at τ V and this case is similar to BD one. (4.8) • Type 2 ( Figure 5(c)): If τ 2 < τ 1 < τ V < τ , the non-BD mode is at τ 2 which is the earliest time. (4.9) We can see that the structure of these equations are the same as the previous ones. The only difference is the upper and lower limit of time integration. It is not hard to generalize to the case where there are L 1 loops later than τ V and L 1 loops earlier than τ V as shown in Figure 4(d).
Although there are two types of diagrams, actually they are unified with the same structure. The diagram in Figure 4 The coefficient is the product of the previous BD one and those contributed by L 1 loops between time τ 0 and τ V , So, the final loop corrected non-BD three-point correlation function is with the loop correction given by It is very interesting to notice that , (4.14) where "1" essentially is just the effective BD coefficients to the lowest order (2.13). So, the ratio between non-BD and BD effects in different sectors are the same. The relevant integrals of 3pt function are: for BD one and for non-BD one, where B ≡ 3λ 2 3200π 3 H 2 . There are exponential oscillation terms in the integrals. In order to find a characteristic scale of initial time, we switch to consider the simpler case by neglecting the exponential oscillations and study the following two integrals , (4.17) .

(4.18)
As for function Q, when τ 0 → −∞ the function in the integral is highly suppressed by the exponential factor. It makes no sense to choose the infinitely past initial time because Q will vanish in that case. Instead, we try to find the conditions for maximal 3pt function. We expect that there exists one initial time τ 0m which can maximize the integral Q.
So far our discussions are based on sub-horizon limit approximations |kτ |, |kτ 0 | 1, but mathematically, it is still meaningful to set τ = 0. This can be justified from the fact that the function Q is very insensitive to the final time τ if |τ 0 | |τ | as well as from Figure 6 where two curves with different initial time are compared. For concreteness, we choose k 1 = k 2 = k 3 = k and τ = 0, then we get the function of τ 0 We need to find dQ dτ 0 τ 0 =τ 0m = 0. Numerically, we can find that the integral is maximized This roughly calibrates the initial time for maximal non-Gaussianity and its corresponding amplitude. The weak dependence on the final time τ and different shape (k 2 /k 1 , k 3 /k 1 ) is shown in the Figure 6. While for function R, it will saturate for early enough initial time (see Figure 6).

Standard result on non-Gaussianity
Before showing the non-Gaussianity under loop corrections, we first review the standard non-Gaussianity, i.e. the tree level result. In the standard procedure, the initial time is chosen to be past infinity. This will cause divergence of the integral. In order to regulate the divergence, we need to adopt the so-called iε prescription which is well understood in standard quantum field theory, but a little problematic in cosmology.
In standard QFT, iε prescription is valid and vital from both mathematics and physics. It not only cures the divergent problems in mathematics, but also ensures that physically the quantum system can evolve from the non-interacting vacuum state in the infinitely past to the true vacuum state with interaction at present.
However, in cosmology, our universe may start from a finite initial time. What's more, de-Sitter inflation has some different non-trivial properties compared to flat space case. The most well-known and serious problem for this iε prescription is the folded limit divergence of non-BD non-Gaussianity. Even for the BD case, iε prescription comes from the scattering problem where particles are initially far away from each other. This is not the case for inflation because the inflation originates from the sub-horizon scales.
Let us first present the result of standard non-Gaussianity. The BD non-Gaussianity can be obtained from the tree level 3pt function calculations by setting τ = 0: where prime denotes that (2π) 3 δ(k 1 +k 2 +k 3 ) is stripped, while the non-BD non-Gaussianity is given by Note that for external line, we only consider the BD modes. Non-BD modes in the external lines essentially are just effectively renormalizing the BD part by contributing term So, we are not going to consider them anymore and alway assume the BD external modes.
The non-Gaussianity can be characterized by the non-Gaussianity shape function F which is defined as where the scale invariance of the correlation functions has been used to show that F only depends on the ratios of different momenta. So, for the standard non-Guassianty, the shape functions for BD and non-BD are given by The corresponding shapes can be seen from Figure 7. Evidently, non-Gaussianity diverges in the folded limit. But this divergence is unphysical. As we will show later, the decay of non-BD modes cures the divergence.

Non-Gaussianity with loop correction
Next, we need to consider the loop corrections. The non-BD non-Gaussianity is maximized roughly when non-BD modes are generated or excited at τ 0m instead of the infinitely past.    Figure 8. BD non-Gaussianity shape function at tree level (left) and loop level (right) for initial time τ 0 = eτ 0m . At tree level, the sharp initial time cut-off will give rise to a fast oscillating non-Gaussianity shape. After including loop correction, the oscillating behavior is suppressed to nearly vanishing value and we nearly recover the usual BD non-Gaussianity shape.
Previously, we only handle the effective value under the sub-horizon approximations, i.e. |kτ | 1. But physically we expect that in the super-horizon case, the decay is quite slow due to the frozen of modes. So, super-horizon and sub-horizon admit completely different behaviors. Tentatively, we can find an intermediate time τ int (|kτ int | 1) to connect these two pieces. When |τ V |> |τ int |, previous sub-horizon approximated results are reliable. While for |τ V |< |τ int |, we can just simply ignore possible loop corrections (which are expected to be very small due to limited time integrations as well as nearly frozen super-horizon modes) and only consider the tree level results with different initial conditions-the renormalized non-BD coefficients at τ int or more explicitly, c k → c eff k (τ int ). So, the loop corrections can be evaluated in the following way where · · · denotes the tree level relevant terms. It is very interesting to note that c eff Non-BD (τ V ) for |τ 0 | |τ |, |τ V |. This suggests that mathematically, the formula for Z loop

Non-BD
can also be used in the super-horizon limit due to its similar behavior. Similar consideration also holds for the BD non-Gaussianities. Based on these arguments, we can still use our previous result derived in the sub-horizon limit to calculate the observable super-horizon non-Gaussianities simply by setting τ = 0: . They can be regarded as the loop corrections or the renormalization factors to the tree diagram. Especially, note that they don't show decay behavior when τ V is pretty small, consistent with our previous physical picture for superhorizon modes.
As we stated before, the loop corrections, in principle, also contain a phase factor exp(iγ) with γ ∼ Bk 5 τ 5 V . We don't consider them because the exact expression is unknown and may be very complicated. The above rough form can be understood from the DRG method. What we want to emphasize is that the amplitude decay is sufficient to suppress the divergence and the fast oscillating phase factor can only be more beneficial due to the dramatic cancellations between positive and negative parts.
Finally, we obtain the BD and non-BD non-Gaussianity shape functions with loop corrections: Next, we give some plots for non-Gaussianity. Recall that in general single field inflation, the power spectrum P ζ = H 2 8π 2 and Σ = H 2 (c s = 1 in our model). The non-Gaussianity estimator f NL = − 10 81 λ Σ = − 10 81 λ H 2 [13]. So, we can express the exponential decay factor in terms of observable quantities as (4.28) We choose parameters P ζ = 10 −9 , f NL = 1 [17]. For non-BD parameters, we use c k = 0.1, θ k = 0. The shape for BD and non-BD non-Gaussianity are shown in Fig. 8 and Fig. 9, Fig. 10. BD non-Gaussianity shape: The non-Gaussianity shape for BD part is shown in Fig. 8. For BD non-Gaussianity, if we choose one initial time sharp cut-off, at tree level the non-Gaussianity shape function shows oscillating behavior due to the oscillating term in the integral. But at loop level, as long as the initial time is not too late which is always the case because the BD starts from very very early time and in principle from nearly past infinity, the oscillating behavior disappears and we nearly recover the usual BD non-Gaussianity shape. Non-BD non-Gaussianity shape: For non-BD non-Gaussianity, it peaks at one specific initial time roughly. Earlier or later initial time can only generate smaller observational non-Gaussianity. What's more, loop corrections cure the folded divergence behavior of non-Gaussianity as we emphasize before. If the initial time is much earlier, which means substantial time for non-BD state to decay, it may be very difficult to observe the remnants of the non-BD information experimentally.
Due to the highly sensitive dependence on initial time, the final shape of non-Gaussianity may show some oscillating features which is not generic and depends on the details of c k , θ k , τ 0 (k) and so on. We expect that in reality, these highly sensitive dependence is fragile and will be averaged or smoothed due to complicated behavior of these functions. The important and generic part is the non-oscillating part with relatively weak dependence on initial time. Note that for BD part, we do not need to use this smoothing functions. Because, in principle, the BD exists from the very early beginning, almost infinitely past. And our exponential correction term is enough to suppress the oscillating parts to get the standard BD non-Gaussianity shape. Nevertheless, for completeness, we provide a typical plot for the non-BD non-Gaussianity shape without smearing initial time.
Therefore, we choose to filter the slowly varying non-oscillating parts by averaging the initial time with a Gaussian distribution centered at τ 0c with width τ 0w : With this smoothing function, the observable non-Gaussianity shape function is Non-BD (4.30) Under Gaussian smoothing, the exponential function e i(k−k 0 )τ 0 will be transformed into a smooth and non-oscillating Guassian function of k centering at k 0 with width 1/τ 0w .
The non-Gaussianity shape of non-BD part (see Fig. 10) includes two features: folded shape peak and squeezed limit shape. The folded shape is mainly contributed by the c k 1 : when k 2 + k 3 − k 1 ∼ 0 or |(k 2 + k 3 − k 1 )τ 0 | π, the sine function in the integral will contribute coherently with oscillations. If we only consider the tree level result, the folded shape value will blow up for very early initial time. But if we include the loop corrected exponential term, its value will be suppressed to nearly vanishing. While the off-diagonal corner (k 2 ∼ 0 part and k 3 ∼ 0 part) shape is contributed by the second and third term in the non-BD non-Gaussianity shape function, namely c k 2 , c k 3 term. When k 2 ∼ 0 and thus k 1 + k 2 − k 3 ∼ k 2 ∼ 0, the sine function c k 3 sin((k 1 + k 2 − k 3 )τ V ) in the integral will contribute coherently as long as |(k 1 + k 2 − k3)τ 0 | π and thus give rise to large value. The arguments are similar to the folded shape one. Similar arguments apply for k 3 ∼ 0.
At tree level, the earlier the initial time, the larger the folded shape peak value. While at loop level, the folded peak will be suppressed to very small value.

Non-interacting limit
In the following part, we are going to consider the non-interacting limit λ → 0 and show that standard results for BD non-Gaussianity can be recovered.
For BD non-Gaussianity, The integral can be evaluated in the following way through integration by parts In the non-interacting limit, D → 0, x 0 → ∞, the above integral gives rise to a factor −2. Thus, we do recover the standard BD non-Gaussianity. Previously, the standard non-Gaussianity is obtained by iε prescription. Here the interaction coupling constant plays the role of ε and regulates the divergence problem. In this sense, we give a natural explanation to the problem of introducing iε in cosmology.
For non-BD non-Guassianity, things are a little more complicated due to the interplay of folded limit and non-interacting limit. Physically, the correct order of taking limits should be like this: fix the coupling strength B first and then examine the non-Gaussianity at folded limit, finally turn off the interactions gradually. After fixing the coupling strength, the initial time is roughly given by τ 0m instead of −∞ in the iε prescription, which can only give a trivial vanishing result in our case. Near the folded limit, we can perform Taylor expansion for δk, then this integral is a normal one and vanishes when δk → 0. So, for a fix coupling strength, the exact folded limit non-Gaussianity vanishes. Then, we take the non-interacting limit, which by continuity also gives rise to a vanishing folded limit non-Guassianity. However for a fixed given coupling strength, globally the amplitude will increase if we let coupling strength goes to zero and choose proper initial time. More specifically, if the shape is not too folded (B(k/δk) 5 is still very small with δk = k 2 + k 3 − k 1 or other permutations), the method for BD one can be applied here. Then, we recover the 1/δk 3 factor in F Non-BD . While the other integral factor (which is −2 for BD one) is highly sentive to the intial time. The scale of amplitude is roughly given by Eq. (4.20), scaling like B −3/5 . This amplitude decreases very quickly once a different initial time is chosen and will vanish if we start from past infinity.
In particular, if we tune the coupling constant λ to be pretty small, for example equivalently let f NL = 10 −6 which is mathematically meaningful nevertheless, we can get non-BD non-Gaussianity which is nearly divergent at the folded limit provided that the proper initial time is chosen. Actually, the amplitude of the folded limit non-Gaussianity is proportional to f −6/5 NL (4.20). However, once the initial time is slightly different, the amplitude decreases dramatically.
The conclusion is that the loop corrections is vital for non-BD three point functions. Even in the extremely weak interacting limit and we can fine tune the initial time delicately to get the ordinary divergent folded non-Guassianity, this apparent divergence breaks down once a little bit different initial time are considered. This is not the feature of previous non-BD three point functions in literature. From this aspect, it is very challenging to observe the imprints of initial non-BD state at present, especially for large k modes.

Conclusion
In this paper, we develop the techniques of calculating one loop diagrams. And we discover a recursion relation which enables us to deal with infinite loop calculations and do resummations. By using these techniques, we show that the decay of non-BD coefficients are consistent with the previous cut-off result except an order one decay factor difference. Our method is enlightening and may shed light on the future loop calculations and resummations in the general context. Furthermore, we analyze the non-Gaussianity under loop corrections. As we expect from the decay of effective non-BD coefficients, the usual divergent non-BD non-Gaussianity at folded limit gets smoothed. What's more, the loop corrected non-BD non-Gaussianities peak at specific initial time and are very sensitive to these initial time. Once we deviate a little bit, these non-Gaussianities will decrease dramatically. So, we conclude that the non-BD non-Gaussianities are very fragile to loop interactions and initial time. Thus, as long as the non-BD state is set up at early enough time, the imprints of these non-Gaussianties on observations may be difficult. These are very different from the previous results in literature where folded limit non-Gaussianities are dominated by the non-BD one due to the divergent behavior.
Besides, we also show that even for BD non-Gaussianity, the loop corrections can have significant influence, playing the role of infinitesimal regulator like iε prescription. The loop correction can not only regulate the divergence problem but also recover the usual result in literature. Thus, loop corrections provide a natural way of introducing iε in cosmology in a more natural and physical way.
Our results are derived based on sub-horizon limit approximations. It is well known that the sub-horizon modes do not feel the presence of gravity much and their behaviour resemble the flat Minkowski space case. So are the loop corrections. One may wonder whether the same properties that we have discussed already exist in flat space quantum field theory. The answer is yes or no. On the one hand, the UV limit of the cosmological perturbations indeed return to flat space quantum field theory. But on the other hand, in usual treatment of flat space quantum field theory, we are interested in the in-out amplitude. Interactions are shut off at the asymptotic past and future. However, in our case, the fluctuations keep on interacting in the asymptotic past. Also, the expansion of the universe exposes anything odd in the UV, if not diluted by the expansion of the universe, to observables at macroscopic scales. Finally, inflation needs a start and may have features, so asymptotic Lorentz symmetry or de Sitter symmetry may not help to determine the vacuum. Those reasons explain the difference between our work and a conventional treatment of flat space quantum field theory.
Although we only consider the general single field inflation and rely on some approximations, the conclusions are expected to apply in more general cases.

B Elliptical coordinate system
Consider two fixed points A, B separated by distance R, and P is the moving point. The distance between A, B, O (the origin) and P is r A , r B , r. Also, assume the angle between OA and OP is θ. The Jacobi matrix is In 3D, the volume element is dV = r 2 sin θdrdθdφ = −rdrd(r cos θ)dφ = −2πrdrdz → 2πr R 3 8 In general, in n-dimensional space, the volume element (after integrating out the angular part) is dV = 2π( √ π) n−3 Γ( n−1 2 ) among them to ensure the regularity of the function. There are 2N coefficients a n , b n , but the regularity near z = 0 will lead to N constraint equations corresponding to each order of Taylor expansion of K in 1/z n . This implies that actually there are only N free parameters. So, we can decompose the above function in terms of the basis functions T n n a n e izu − b n z n = n A n T n (z) , (C.6) We can express A n in terms of a n , b n by solving the above equation. Also, note that in the basis functions T n , 1 z only appears in T 1 , so we can get a very important relation