Strong restriction on inflationary vacua from the local gauge invariance III: Infrared regularity of graviton loops

It has been claimed that the super Hubble modes of the graviton generated during inflation can make loop corrections diverge. Even if we introduce an infrared (IR) cutoff at a comoving scale as an ad hoc but a practical way for the regularization, we encounter the secular growth, which may lead to the breakdown of perturbative expansion for a sufficiently long lasting inflation. In this paper, we show that the IR pathology concerning the graviton can be attributed to the presence of residual gauge degrees of freedom in the local observable universe as in the case of the adiabatic curvature perturbation. We will show that choosing the Euclidean vacuum as the initial state ensures the invariance under the above-mentioned residual gauge transformations. We will also show that as long as we consider a gauge invariant quantity in the local universe, we encounter neither the IR divergence nor the secular growth. The argument in this paper applies to general single field models of inflation up to a sufficiently high order in perturbation.


I. INTRODUCTION
The inflationary spacetime leads to the generation of gravitational waves. Even though the amplitude of gravitational waves is smaller than the amplitude of the adiabatic curvature perturbation, the detection of the primordial gravitational waves generated during inflation is expected to be within our reach. The measurement of the primordial gravitational waves is crucially important to uncover the model of inflation, providing information which cannot be obtained through the measurement of the scalar perturbations. In particular, the amplitude of the gravitational waves can directly measure the energy scale of inflation unlike the amplitude of the adiabatic curvature perturbation, which is sensitive also to the detailed dynamics of the inflationary universe. We can find an example which highlights the importance of measuring the primordial gravitational waves also for non-linear perturbations, say in Refs. [1][2][3], which elucidated the impact of parity violation in gravity sector on the bi-spectrum of the primordial gravitational waves. Thus, it will be profitable to study the method to predict the primordial gravitational waves also in the presence of the non-linear interaction. In this paper, our focus is on the loop correction due to the primordial gravitational waves.
It is known that loop corrections of a massless scalar field generated during inflation can suffer infrared (IR) divergences . Since a massless scalar field yields the scale-invariant spectrum in the IR limit as P(k) ∝ 1/k 3 , a naive loop integral yields a factor d 3 k/k 3 ∼ dk/k, which logarithmically diverges. As is expected from the fact that the mode equation of the tensor perturbation, which we refer to also as the graviton, is nothing but the mode equation for a massless scalar field, graviton loop corrections also appear to yield IR divergences. To quantify the graviton loop corrections, we need to provide a way to regularize them. One may propose to introduce an IR cutoff say at a comoving scale k IR as a practical way for the regularization. However, this will not provide a satisfactory solution, because the loop integral of the super Hubble modes gives aH kIR dk/k ∼ ln(aH/k IR ), which logarithmically increases in time. Here a and H are the scale factor and the Hubble parameter of the background spacetime, respectively. Compared with the tree level contribution, the loop corrections are typically suppressed by the Planck scale as (H/M pl ) 2 with M 2 pl ≡ (8πG) −1 . However, this suppression may be compensated by the secular growth for a sufficiently long lasting inflation. (A thorough overview about the possible origins of the IR divergences can be found in the review paper [43] by Seery.) The IR behaviour of the graviton has attracted attention also as a possible origin of the running of coupling constants [10-12, 24, 44-53]. Tsamis and Woodard claimed that the logarithmically growing secular effect due to the graviton loops can lead to the screening of the cosmological constant, suggesting that the screening may provide a dynamical solution to the cosmological constant problem [24]. More recently, Kitamoto and Kitazawa studied the IR effect on the gauge coupling and claimed that the secular growth from the graviton loops can screen the gauge coupling [34,35]. A related issue is discussed for the U(1) gauge field in Refs. [36,37]. If the secular growth due to the graviton loops is actually physical, it will provide a phenomenological impact. However, these secular growths originate from the increasing IR contributions, and hence the predicted secular effects significantly depend on the regularization method of IR contributions. Therefore, to trust the secular growth as a physical effect, its presence should be verified based on a rigorous method for the regularization. The graviton loop corrections have been mostly discussed in the exact de Sitter background and their regularity is still under debate even at the linear level. In Ref. [54], Higuchi et al. claimed the existence of a regular two-point function (see also Refs. [55,56]), while Miao et al. objected against it in Ref. [57]. This issue was discussed also in Refs. [58][59][60].
The IR pathology has been studied more intensively for the adiabatic curvature perturbation [61][62][63][64][65][66][67][68][69][70][71][72][73][74][75][76]. In the presence of gravity, we need to remove the influence of gauge degrees of freedom in calculating the observable fluctuations. In the cosmological perturbation theory, gauge artifacts are usually removed by employing gauge conditions that completely fix the coordinate choice. However, when we consider actual observations, the gauge conditions may not be fixed in a suitable way such that preserves the genuine gauge invariance. In fact, we can observe the fluctuations only within the region causally connected to us, which is limited to a finite portion of each time constant slice. Therefore, precisely speaking, the gauge conditions used in the conventional cosmological perturbation theory, which require the information all over the time constant slice of the universe, does not match the actual process of observations. To regularize the IR contributions for the curvature perturbation, it is necessary to take into account this subtle issue [62,63,[65][66][67]. Since in actual observations we can impose the gauge conditions only in the observable region, there inevitably appears the ambiguity associated with the degrees of freedom in choosing coordinates in the outside of the observable region. Such residual coordinate degrees of freedom can be attributed to the degrees of the freedom in the boundary conditions of our observable local universe. It was shown that requesting the invariance under the change of such residual coordinate degrees of freedom in the local universe can ensure the IR regularity and the absence of the secular growth [62,63,[65][66][67]. This is, so to speak, because we can absorb the singular IR contributions by the residual coordinate degrees of freedom (see also Refs. [77,78]). This issue is reviewed in Ref. [67]. It has been pointed out that the residual coordinate degrees of freedom also can affect the IR behaviour of the graviton [63,67]. In this paper, we will show that, when we require the invariance under the residual coordinate transformations, the IR regularity and the absence of the secular growth are guaranteed also for the graviton loops. (The IR issues about the entropy perturbation were studied in Refs. [79,80] and those about a test scalar field in a de Sitter background were studied in Refs. [81][82][83][84][85][86]. ) For the graviton, the relation between the IR divergence and the gauge artifact has been discussed during the past few decades. Allen showed that the IR divergence in the free graviton propagator which appears for some particular values of the gauge parameter can be understood from the fact that the gauge fixing term does not select a unique gauge for these particular values of the gauge parameter [87]. (See also Refs. [88][89][90][91].) Even if we properly choose the gauge parameter, it is known that the transverse traceless mode can still suffer from the IR divergence through the loop corrections, which is the target of this paper. The connection between the IR divergence and the gauge artifact was pointed out also in Ref. [92], where the secular growth predicted by Tsamis and Woodard in Ref. [24] was reexamined. It was shown that the spatially averaged Hubble expansion computed by Tsamis and Woodard is not invariant under the change of the time slice and hence the screening effect which shows up in their averaged Hubble parameter suffers from the gauge artifact [92,93]. Meanwhile, in Ref. [92], it was shown that the locally defined Hubble expansion which may mimic the observable Hubble rate becomes time-independent. This example suggests that computing observable quantities, unaffected by the influence from the outside of the observable region, will play an important role to quantify the IR effects of the graviton [92] (see also Ref. [94]). This paper is organized as follows. In Sec. II, we will briefly show how the IR divergence and the secular growth can appear when we naively compute the loop corrections for the curvature and graviton perturbations. Then, in Sec. II B, we will point out the presence of the residual coordinate degrees of freedom in the local observable universe, which describes the influence from the outside of the observable region. One way to preserve the invariance under the residual coordinate transformations will be discussed in Sec. II C. Using the prescription which will be introduced in Sec. III, we will show, in Sec. IV and Sec. V, that when we choose the Euclidean vacuum as the initial state, the IR regularity and the absence of the secular growth are ensured. In this paper, we will discuss the inflationary universe which contains a single scalar field and we will not directly discuss the pure gravity setup, although our argument may also provide some insight into the latter case.

II. OVERVIEW OF IR ISSUES
In this section, we will overview the IR issues of the curvature perturbation and the graviton perturbation, which are the main target of this paper. In this paper, we consider a standard single field inflation model whose action is given by where M pl is the Planck mass and φ is the dimensionless scalar field which is obtained by dividing an ordinary scalar field by M pl . We choose the time slicing by adopting the uniform field gauge Under the ADM metric decomposition, which is given by we further decompose the spatial metric g ij as where a ≡ e ρ is the scale factor, ζ is the so-called curvature perturbation and δγ ij is the graviton perturbation which satisfies the transverse and traceless condition Since δγ ij is traceless, we find det γ = 1.

A. Various divergences from the curvature and graviton perturbations
In this subsection, after we briefly review the linear perturbation theory for the scalar and tensor perturbations, we will summarize several different origins of the divergences that possibly appear in the loop corrections of these two perturbations.

Scalar perturbation
The quadratic action for the scalar perturbation ζ, which describes the evolution of the interaction picture field ζ I , is given by and its equation of motion is given by Here, dot denotes the differentiation with respect to cosmological time t. Here, for notational convenience, we introduce the horizon flow functions as with n ≥ 2. However, we do not assume that these functions are small, and leave the background inflation model unrestricted.
For quantization, we expand the curvature perturbation ζ I (x) as where a k is the annihilation operator, which satisfies The mode function v s k (t) should satisfy 11) and is normalized as where the Klein-Gordon inner product is defined by Notice that Eq. (2.11) states that v s k (t) becomes time-independent in the IR limit k/(e ρρ ) ≪ 1. When we choose the mode function for the adiabatic vacuum, which approaches the WKB solution in the UV limit k/(e ρρ ) ≫ 1, the power spectrum becomes almost scale invariant in the IR limit as (2.14) where we evaluated v s k (t) at the Hubble crossing time t = t k with k = e ρ(t k )ρ (t k ), since the curvature perturbation gets frozen rapidly after t k .
When we assume that the corresponding free theory has an almost scale-invariant spectrum in the IR limit, a naive consideration can easily lead to the IR divergence due to loop corrections. For instance, one may expect that an interaction vertex which contains the curvature perturbation without the derivative operator yields the factor whose momentum integral logarithmically diverges in the IR as d 3 k/k 3 for the scale-invariant spectrum. Even if the spectrum is not exactly scale invariant as in Eq. (2.14), the deep IR modes give a significant contribution to Eq. (2.15). Following Ref. [67], we refer to such unsuppressed contribution due to the integral over small k as the IR divergence (IRdiv). When we introduce an IR cutoff for the regularization, say, at the Hubble scale at the initial time t i , the super Hubble (superH) modes in the variance of ζ I give rise to the secular growth which is logarithmic in the scale factor a = e ρ as Then, the loop corrections, which are suppressed by an extra power of the amplitude of the power spectrum (ρ/M pl ) 2 , may dominate in case the inflationary epoch lasts sufficiently long, leading to the breakdown of the perturbative expansion. We refer to the modes with e ρ(ti)ρ (t i ) k e ρ(t)ρ (t) as the transient IR (tIR) modes and refer to the enhancement of the loop contributions due to the tIR modes as the IR secular growth (IRsec). From the definition, it is clear that the tIR modes were in the sub Hubble (subH) range at the initial time t i , but have been transmitted into the superH ones by the time t. The influences of the IRsec have been discussed by introducing an IR cutoff in Refs. [95][96][97][98][99][100][101][102]. We refer to the case when both the IRdiv and the IRsec are absent as IR regular.
So far, we discussed the secular growth which originates from the momentum integration with the time coordinates of the interaction vertices fixed. However, the time integration also can yield the secular growth. If the contribution from the interaction vertex in the far past remains unsuppressed, it will diverge when we send the initial time to the infinite past. We refer to the secular growth due to the temporal integral as the SG, discriminating it from the previously discussed IRsec. (Regarding the SG, see also Refs. [103,104].)

Graviton perturbation
The quadratic action for the graviton perturbation δγ ij , which describes the evolution of the interaction picture field δγ ij I , is given by and its equation of motion is given by We quantize δγ ij I as where λ labels the helicity, e ij (k, λ) is the transverse and traceless polarization tensor which satisfies and a (λ) k is the annihilation operator which satisfies When the graviton perturbation is isotropic, its variance (in the coincidence limit) is given by where P ij is the transverse traceless projection tensor: and P t (k) is the power spectrum of the graviton perturbation: Here, the factor 2 counts the number of helicity. Since the equation for v t(λ) k is identical to the one for a massless scalar field, the graviton spectrum in the adiabatic vacuum is almost scale invariant in the IR limit as Integrating over the angular coordinates in Eq. (2.22), we obtain Similarly to the curvature perturbation, we find that the IR and tIR modes in Eq. (2.26) yield the IRdiv and IRsec, respectively. The influence of the IRsec due to graviton loops is discussed in Ref. [98]. Meanwhile, the interaction vertices with the graviton also can yield the SG in the same way as the curvature perturbation.

Solving the constraint equations
Eliminating the Lagrange multipliers N and N i , we derive the action written solely in terms of the dynamical fields ζ and δγ ij [105]. In the gauge defined by Eqs. (2.2) and (2.5), the constraint equations are given by where D i is the covariant differentiation associated with the spatial metric g ij , and are the extrinsic curvature and its trace, respectively. We expand the metric perturbations as where we use the subscript I to express the first order curvature and graviton perturbations, since they correspond to the interaction picture fields. Then, the n-th order Hamiltonian constraint and the momentum constraints are expressed in the form where ∂ 2 denotes the Laplacian. For n = 1, H 1 and M i 1 are 0 and for n ≥ 2, H n and M i n consist of n interaction picture fields (either ζ I or δγ ij I ).
Since the constraint equations (2.34) and (2.35) are elliptic-type equations, we need to employ a (spatial) boundary condition to fix a solution of N n and N i n . A general solution of (2.34) and (2.35) in the absence of the graviton perturbation can be found in Appendix of Ref. [66]. An extension to include the graviton perturbation proceeds in a straightforward manner and the general solution is given as where G n (x) and G i n (x) are arbitrary solutions of the Laplace equations and n-th order in perturbation. Since the function G i n (x) contributes only through its transverse part, the number of introduced independent functions is three. By employing appropriate boundary conditions at the spatial infinity, the boundary conditions for these elliptic-type equations will be uniquely fixed. Substituting thus obtained expression for N and N i , the action S = d 4 xL[ζ, δγ ij , N, N i ] can be, in principle, expressed only in terms of the dynamical fields ζ and δγ ij . Then, the evolution of ζ and δγ ij is governed by the non-local action with the inverse Laplacian.
As was pointed out in Refs. [62,63], the inverse Laplacian ∂ −2 may enhance the singular behavior of perturbation in the IR limit by introducing a term with the factor 1/k 2 where k is a comoving momentum of the constituent fields. (A detailed explanation is given in the review article [67].) We will return to this issue in Sec. IV B 2.

Observable local patch
To discuss the observable quantities, we introduce the observable region as the region causally connected to us. We express the observable region on the time slice at the end of inflation t f and its comoving radius as O t f and L t f , respectively. The causality requires that L t f should satisfy where t 0 is the present time. The cosmological observations can measure the n-point functions of the fluctuation with the arguments (t f , x) contained within the observable region O t f . For later use, we refer to the causal past of O t f as the observable region O and refer to the intersection between O and a t-constant slice Σ t as O t . We approximate the comoving radius of the region O t as .
As L t is approximated by L t ≃ 1/(e ρρ ) in the distant past, the effects of the superH modes with k e ρ(t)ρ (t) can be understood as the influence from the outside of the observable region O. These modes potentially affect the fluctuations in O t f by two means. One is due to the non-local interaction through the inverse Laplacian ∂ −2 , while the other is through the Wightman functions and Even if the spatial distance |x 1 − x 2 | is bounded from above by confining both x 1 and x 2 within the observable region, the contribution to the Wightman functions from the superH modes with k ≤ |x 1 − x 2 | −1 is not suppressed. The superH modes make these Wightman functions divergent for scale-invariant or red-tilted spectrum. (See Eqs. (2.16) and (2.26).) The regularity of the contribution from the superH modes can be verified, only if their contribution is suppressed by an additional factor of k.

Residual gauge degrees of freedom
In Sec. II B 2, we introduced the observable region O, which is a limited portion of the whole universe. We claim that the observable quantities must be composed of the fluctuations within O. Since only the fluctuations within O are relevant, there is no reason to request the regularity of the fluctuations at the spatial infinity in solving the elliptic-type constraint equations (2.34) and (2.35) (at least, at the level of Heisenberg equations of motion). Then, there arise degrees of freedom in choosing the boundary conditions, which are described by arbitrary homogeneous solutions of the Laplace equation, G n (x) and G i n (x) in Eqs. (2.36) and (2.37). These arbitrary functions in N and N i can be understood as the degrees of freedom in choosing the coordinates. Since the time slicing is fixed by the gauge condition (2.2), the residual gauge degrees of freedom can reside only in the spatial coordinates x i .
Let us consider these residual coordinate transformations associated with G n (x) and G i n (x). We add the subscript gl to the original global coordinates, in order to reserve the simple notation x for the coordinates after transformation. As we have shown in Refs. [62,63], the coordinate transformations xgl → x are specified as We should note that, once we substitute the expressions for N and N i to obtain the equation of motion solely written in terms of the curvature perturbation ζ and the graviton perturbation δγ ij , the symmetry under the residual coordinate transformations is lost, because N and N i depend on the specified boundary conditions. In this sense the coordinate transformations (2.42) should be distinguished from the standard gauge transformations that leave the overall action invariant. Therefore, to avoid confusion, we distinguish the coordinate transformations (2.42), writing it in the italic font as the gauge transformation.
Among the residual gauge transformations, we focus on the one with where s(t) is a time-dependent function and S ij (t) is a time-dependent traceless tensor. When we consider the time-dependent dilatation transformation, the homogeneous part of the curvature perturbation transforms as ζ → ζ − s(t) .
(2.44) (Precise meaning of this transformation will be explained later.) In Refs. [65,66], we showed that preserving the invariance under the dilatation transformation parametrized by s(t) is crucial to show the regularity of the loops of the curvature perturbation. Intriguingly, the transformation (2.43) shifts the graviton perturbation as at the linear level, which is analogous to the shift for ζ. Although the non-linear extension of the above transformation is rather non-trivial, this observation suggests that analogous proof of the IR regularity may also work for graviton loops. The relation between the IRdiv due to graviton loops and the homogeneous shift (2.45) was pointed out several times. Gerstenlauer et al. [69] and Giddings and Sloth [70] showed that the leading IRdiv of the graviton loops can be attributed to the change of the spatial coordinates in the form (2.43) with s(t) = 0 due to the accumulated effect of the IR graviton.

C. Genuine gauge invariance and quantization
The observable fluctuations should not be affected by the residual gauge degrees of freedom, which were discussed in the preceding subsection. In this subsection we discuss how to introduce a quantity which is invariant under the residual gauge transformations. We call such a quantity a genuinely gauge invariant quantity. One may think that the genuine gauge invariance will be preserved by fixing the residual gauge degrees of freedom completely. If we could perform a complete gauge fixing by employing appropriate boundary conditions for N and N i at the boundary of the observable region O, the IRdiv and IRsec will not appear, because the maximum wavelengths of fluctuations in such a gauge would be bounded approximately by the size of O. We pursued this possibility in Refs. [61,67]. When we perform the quantization after complete gauge fixing or equivalently perform the quantization within the local observable region O, the global translation invariance is not easily preserved any longer, leading to technical complexities. To avoid the complexities, in Ref. [61], first we performed the quantization in the whole universe, and then we fixed the coordinates by carrying out the residual gauge transformation. In this manner we showed that the absence of the IRdiv and IRsec is guaranteed if the initial fluctuation does not suffer from these IR pathologies. However, it turned out that the IRdiv can arise even in the initial fluctuation after we perform the residual gauge transformation to employ complete gauge fixing [67].
Here, following Refs. [62,63], we perform the quantization, taking into account the whole universe without fixing the residual gauge degrees of freedom, which allows us to keep the global translation invariance manifestly. Then, we construct a genuine gauge invariant operator and choose an initial state which will be understood as the genuine gauge invariant state. Since the time slice is uniquely specified by the gauge condition (2.2), the genuine gauge invariance will be ensured when a quantity preserves the invariance under a change of spatial coordinates.
To construct a genuinely gauge invariant operator, we consider correlation functions for the scalar curvature of the induced metric on a φ =constant hypersurface, s R. Since s R itself transforms as a scalar quantity under spatial coordinate transformations, the correlation functions of s R are not invariant. However, the n-point function of s R will become gauge invariant, if we specify its n arguments in a coordinate-independent manner. The distances measured by the spatial geodesics that connect all the pairs of n points characterize the configuration in a coordinate independent manner. Here, we adopt a slightly easier way, specifying the n spatial points by the geodesic distances and the directional cosines that are measured from an arbitrarily chosen reference point. Although the choice of the reference point and the frame is a part of the residual gauge , this ambiguity will not matter as we will choose a quantum state such that does not break the spatial homogeneity and isotropy of the universe.
Our geodesic normal coordinates are introduced by solving the spatial geodesic equation on each time slice: where s Γ i jk is the Christoffel symbol for the three dimensional spatial metric e 2ζ γ ij , and λ is the affine parameter. The initial "velocity" is given by A point x in the geodesic normal coordinates is identified with the end point of the geodesic x i gl (x, λ = 1) in the original coordinates. Perturbatively expanding x i gl in terms of x i , we obtain Notice that the relation between xgl and x depends on the metric perturbations, which become quantum operators after quantization. Finally, we find that, by means of the geodesic normal coordinates, the genuinely gauge invariant variable is given by In order to calculate the n-point functions of g R, we also need to specify the quantum state such that is invariant under the residual gauge transformations. However, in the present approach, we cannot directly discuss this invariance as a condition for allowed quantum states. This is because the residual gauge degrees of freedom cease to exist, when we quantize fields in the whole universe. Let us recall the discussion in the case of the curvature perturbation ζ [65][66][67]. By construction, the operator g R is not affected by the residual gauge degrees of freedom. However, the n-point functions of g R can be correlated to the fields in the causally disconnected region. In Sec. II B 3, we discussed two ways by which the outside of the observable region O can affect the fluctuation in O. One is through the boundary conditions for the inverse Laplacian ∂ −2 . Since changing these boundary conditions is nothing but performing the residual gauge transformation (see Sec. IV B 2), the n-point functions of the genuine gauge invariant curvature perturbation g R are not affected even if we restrict the integration region of ∂ −2 to the region O. Therefore, as long as we consider a genuinely gauge invariant operator, the inverse Laplacian ∂ −2 never gives the conjunction between the inside and the far outside of O. The other leak of the influence from the outside of O is due to the long-range correlation through the Wightman functions, which can remain even if we consider genuinely gauge invariant variables. Therefore this long-range correlation may give a possible origin of the IRdiv and IRsec. In the case with the curvature perturbation ζ, it is shown that requesting the IR regularity by suppressing the long-range correlation constrains the quantum state of the inflationary universe [65]. Interestingly, the IR regularity condition on quantum states can be interpreted as the condition that requests the quantum states to be unaffected by the time-dependent dilatation transformation, one of the residual gauge degrees of freedom [65]. In Appendix A, we show that also for the graviton perturbation, a similar genuine gauge invariance condition on quantum states is derived from the IR regularity condition.

III. PRELIMINARIES
In this section, as a preparation to analyze the n-point functions of the genuine gauge invariant curvature perturbation, we introduce a family of canonical variables. First, in Sec. III A, we describe the basic formulation for the canonical quantization in terms of the original set of variables ζ, δγ ij and their conjugate momenta. In Sec. III B, we introduce a family of alternative sets of canonical variables, in terms of which the proof of the IR regularity becomes more transparent.

A. Canonical quantization
For notational simplicity, we suppress the subscript "gl" in this subsection. In the following discussion, we express the action for the curvature perturbation ζ and the graviton perturbation δγ ij derived by solving the Hamiltonian and momentum constraint equations as which includes the non-local integration operator ∂ −2 . Here, L dyn denotes the functional form of the Lagrangian density obtained after we eliminate the Lagrange multipliers N and N i . We also introduce the Hamiltonian H and the Hamiltonian density H as where we introduced the conjugate momenta as .
This set of canonical variables Φ ≡ {ζ, π, δγ ij , π ij } should satisfy the standard commutation relations and is the tensorial delta function with the transverse traceless projection.

B. Canonical transformation associated with residual gauge transformations
In this subsection, we introduce a family of alternative sets of canonical variables whose HamiltonianH(t) is written only in terms of where s(t) and S ij (t) are arbitrary time-dependent function and symmetric-traceless matrix, respectively. We treat both s(t) and S ij (t) perturbatively, assuming that they are as small as metric perturbations. We also show that s(t) and S ij (t) without time differentiation are contained in the HamiltonianH(t) only in the combination described in Eq. (3.8).

Introducing new canonical variables
For illustrative purpose, we first consider a coordinate transformation with s and S ij time independent, which induces constant shiftsζ(x) − s and δγ ij (t, x) − S ij . To be concrete, we consider the coordinate transformation xgl → x with and For the time being, we will not specify the terms of O(S 2 ) in Λ i T j . Next, we consider the change of the spatial metric g ij under the gauge transformation (3.9). As is addressed in Ref. [65], when we neglect the graviton perturbation, setting Λ i T j = 0, the dilatation transformation changes the spatial metric as where we have definedζ(x) ≡ ζ(xgl ). We find that the curvature perturbation ζ(xgl ) transforms toζ(x) − s, which suggests that this scaling transformation can be used to find the canonical variablesΦ that are subjected to the necessary constant shift. Compared with the curvature perturbation, finding a transformation which shifts the graviton perturbation by −S ij is much more non-trivial, particularly at non-linear order. Therefore, introducing the transverse traceless tensor δγ ij , we express the spatial metric obtained after the coordinates transformation (3.9) as with the requested shift by −S ij . In the following, we assume that s(t) and S ij (t) are of the same order asζ and δγ ij . From Once the functional form of Λ i T j is determined, Eq. (3.13) specifies δγ ij order by order in perturbation. By expanding the inverse matrix of Λ i T j as (3.14) On the right hand side, we explicitly wrote only the linear order in perturbation. Since the left hand side of Eq. (3.14) is independent of S ij , the field δγ ij should be defined so that the S ij dependence on the right hand side is canceled. In particular, Eq. (3.14) states that δγ ij (xgl ) should agree with δγ ij (xgl ) at the linear order in perturbation. Since the diffeomorphism invariance of the action implies that the Lagrangian densities for g ij andg ij should take the same functional form, using the Lagrangian density L dyn in Eq. (3.1), we can express the action forg ij as Next, we extend the above argument to time-dependent transformations with where Λ i T j (t) is a functional of S ij (t) whose explicit form will be specified later. Similarly to the case of constant Λ i j , we introduce a new set of canonical variablesΦ bỹ with the formal definition of the conjugate momenta given bỹ .
In general, the residual gauge transformation is not well-defined in the whole universe, since the transformation can diverge at the spatial infinity. However, the residual gauge transformations (3.16), exceptionally, keeps the variables finite in the whole universe as is manifest from the above relations. Therefore, we can consistently discuss quantum theory in terms of the set of canonical variablesΦ as well.

Commutation relations
Next, we will show that the variablesΦ = {ζ,π, δγ ij ,π ij }, defined in Eqs. (3.17), (3.18) and (3.19), satisfy the standard commutation relation. Because of the time variation of Λ i j , the partial time derivative with the original global spatial coordinates xgl fixed differs from the one with the new coordinates x fixed. We choose the transformation matrix Λ i T j such that the difference between the two partial time derivative operations does not give S ij (t) without the time derivative. Then, we find 20) or equivalently In fact, with this choice of Λ i T j , we obtain where 1 denotes the unit matrix and the round brackets on indices represent symmetrization. The terms with spatial derivative on the right hand sides of Eqs.  Using Eq. (3.22), we find that the conjugate momentumπ is related to π as In the second equality we used which is derived by changing the spatial coordinates in the action from x to xgl . Deriving the relation betweenπ ij and π ij is more non-trivial, but using where in the second equality we used and which can be derived by using Eq. (3.23), we obtaiñ Here, we simply assume that the operator ordering is properly chosen.
Once we establish the relations between the two sets of the canonical variables, Φ andΦ, the commutation relations for Φ yield the commutation relations forΦ. Using Eqs. (3.4), (3.17) and (3.26), we obtain Similarly, using Eqs. (3.5), (3.18) and (3.31), we find In the second equality we noted the tensorial delta function δ by expanding δγ ij (xgl ) as and we used the factor e −3s(t) to change the argument from xgl to x. The remaining commutation relations can be shown in a similar way and hence we can verify thatΦ is actually qualified as a set of canonical variables.
Solving Eq. (3.20), we can determine the transformation matrix Λ i T j (t). As a boundary condition to solve the first order differential equation, we employ the condition at the end of inflation t f . As we have chosen Λ i T j (t) so as to satisfy det Λ T (t f ) = 1, Eq. (3.24) guarantees that the condition (3.25) holds for all t. Notice that we can formally solve Eq. (3.20) as using the time ordered product denoted by the operator T . Perturbatively expanding Λ i T j (t) with respect to S ij (t) to the next to leading order, we obtain (3.36)

Hamiltonians
Next, we compute the Hamiltonian forΦ defined bỹ where H(x) is the Hamiltonian density defined in Eq. (3.2). Rewriting the graviton part is slightly non-trivial, but this can be confirmed as follows. When we express H(t) in terms of these two variables transform as standard tensors in three dimensions intõ leaving aside the factor e 3s(t) , which will be absorbed to make the combinationζ(x) − s(t) (see Ref. [66]). Then, since γ ij and γ ij are given by [exp(δγ)] ij and [exp(δγ − S)] ij , respectively, we can verify Eq. (3.39). Next, collecting the quadratic terms inΦ from the Hamiltonian, we identify the non-interacting Hamiltonian, which is exactly the same form as the one for Φ. Since bothζ and δγ ij , which are massless fields, appear with spatial derivative operators in the non-interacting Hamiltonian, the shifts by −s(t) and −S ij (t), respectively, are eliminated. We also introduce the interaction Hamiltonian as withH where is the interaction Hamiltonian density for Φ.

Short summary of the strategy
It will be instructive to note the two important properties of the interaction Hamiltonian densityH I , given in Eq. (3.42), which will become crucial in the discussion of the IR regularity: First, the fieldsζ(x) and δγ ij (x) always accompany the time dependent parameters s(t) and S ij (t) asζ(x)−s(t) and δγ ij (x)−S ij (t). Second, s(t) and S ij (t) which are not accompanied bỹ ζ(x) and δγ ij (x), respectively, are always differentiated with respect to time. When we consider only the adiabatic perturbation, we can provide the new set of canonical variables which fulfills these two properties simply by considering the time-dependent dilatation transformation xgl → x with xgl = e −s(t) x, which is one of the residual gauge transformation [65,66]. At the linear order, the residual gauge transformation with Eq. (2.43) shifts the spatially homogeneous part of the graviton perturbation. However, at non-linear order, due to the non-commutativity between matrices, the residual gauge transformation with (2.43) does not immediately introduce the shift of all the graviton perturbation in the action.
To provide the new set of canonical variables with their homogeneous parts shifted, we introduced a more non-trivial transformation (3.16). By choosing δγ ij as in Eq. (3.18), the first property can be ensured. To guarantee the second property, we determined (the time dependence of) Λ i T j (t), requesting Eq. (3.20). Then, s(t) and S ij (t) without the time derivative appear neither on the right hand sides of Eqs. (3.22) and (3.23) nor in the HamiltonianH(t). Using this property, we later show the IR regularity of graviton loops in a parallel way to the case of the curvature perturbation.

C. Coarse grained gauge invariant operator
In Sec. II C, we introduced the genuine gauge invariant variable g R, using the geodesic normal coordinates. Changing the spatial coordinates to the geodesic normal coordinates also modifies the UV contributions. Tsamis and Woodard [106] pointed out that using the geodesic normal coordinates can introduce an additional origin of UV divergence, yielding contributions which may not be renormalized by local counter terms [107]. This is presumably because specifying the spatial distance precisely in the presence of the gravitational perturbation requires taking account of all short wavelength modes. In any realistic observations, what we actually observe is a smeared field with a finite resolution. However, it is not so trivial how to describe a realistic smearing in a genuinely gauge invariant manner. Here, in order to remove the UV contribution in the measurement of the position, we replace the geodesic normal coordinates with approximate ones which are not affected by the UV contributions.
In this paper, we compute the n-point functions at the end of inflation t = t f . Then, in place of the geodesic normal coordinates, we use the "smeared" coordinates x i which are related to the global coordinatesx i gl as 1 where we replaced s(t f ) and S ij (t f ) in the transformation matrix Λ i j (t f ) with the smeared metric perturbations, Here, W t (x) is the window function which takes the non-vanishing value in the local region O t and which implicitly depends on the values of However, Λ i T j at t = t f is exceptionally determined by the value of S ij only at t = t f owing to the boundary condition (3.34). Although gζ and δ gγ ij appear on the right-hand sides of Eqs. (3.45) and (3.46), we can iteratively define gζ and δ gγ ij by these expressions. Using the quantities introduced above, we define g ζ(t f , x) and δ g γ ij (t f , x) as Notice that xgl depends only on ζ and δγ ij but is independent of their canonical conjugate momenta. Hence, we can define gζ (t f ), δ gγ ij (t f ), g ζ(t f , x), and δ g γ ij (t f , x) without ambiguity of the operator ordering. The fields g ζ(t f , x) and δ g γ ij (t f , x) introduced above are not genuinely gauge invariant. However, we can show that R x g ζ(t f , x) and R x δ g γ ij (t f , x) preserve the invariance under the particular residual gauge transformation given in Eq. (3.44), where , · · · (3.50) represents an operator such that manifestly suppresses the IR contributions by acting on the fields g ζ(t, x) and δ g γ ij (t, x). Here, the subscript associated with R x specifies the argument of the fields on which the operator acts. In the following, we will address the IR regularity of the n-point functions of R x g ζ(t f , x) and R x δ g γ ij (t f , x).

IV. EUCLIDEAN VACUUM AND REGULARIZATION SCHEME
In order to compute genuinely gauge invariant quantities, we also need to specify the quantum state so as not to be affected by the residual gauge degrees of freedom. However, as we mentioned in Sec. II C, the genuine gauge invariance of the quantum state cannot be directly discussed in our current approach. Hence, focusing on the invariance under the restricted class of transformations (3.16), we discuss the equivalence among quantum states specified in terms of the set of variablesΦ with various values of s(t) and S ij (t). As is discussed in Ref. [66], the boundary condition of the Euclidean vacuum selects the unique quantum state irrespectively of the choice of canonical variables connected by the dilatation transformation. Namely, as long as we choose the Euclidean vacuum, the quantum state is unaltered by the dilatation scaling. In this section, we extend this argument to include the graviton perturbation, using different sets of the canonical variables Φ andΦ, introduced in the previous section by performing the residual gauge transformation from x i gl to x i . We will show that, employing the boundary condition of the Euclidean vacuum, we can select the unique quantum state irrespective of the choice of the canonical variables. In Sec. IV B, using this property of the Euclidean vacuum, we will reformulate the perturbative expansion to give a proof of the IR regularity.

A. Euclidean vacuum
In this subsection, we briefly summarize the basic properties of the Euclidean vacuum. In the case of a massive scalar field in de Sitter space, the boundary condition specified by rotating the time path on the complex plane can be understood as requesting the regularity of correlation functions on the Euclidean sphere which can be obtained by the analytic continuation from those on de Sitter space. The vacuum state defined in this way is called the Euclidean vacuum state. Here, we denote by the Euclidean vacuum the state specified by a similar boundary condition in more general spacetime.
To be more precise, we define the Euclidean vacuum by requesting the regularity of the n-point functions, where a = 1, · · · , n and T c represents the path ordering along the closed time path, For simplicity, here we assume that e ρ(t)ρ (t) is rapidly increasing in time so that We added the subscript E to the expectation values for the Euclidean vacuum defined in terms of the canonical variables Φ.
For the canonical variablesΦ, the boundary condition of the Euclidean vacuum is similarly given by The Euclidean vacuum is expected to be invariant under the residual gauge transformations, since the above boundary conditions of the Euclidean vacuum are formally independent of the choice of canonical variables. In fact, we can show the equivalence between the expectation values, where the operators O andÕ are related with each other by the relations (3.17), (3.18), and (3.19). A more detailed explanation regarding the uniqueness of the Euclidean vacuum can be found in Ref. [66] and the argument there can be extended to include the graviton modes in a straightforward manner. We will find that the distinctive property (4.5) is crucial in showing the IR regularity for the Euclidean vacuum.

B. Rewriting the n-point functions
In this subsection, we rewrite the expression for the n-point functions into a more suitable form to examine the regularity of the IR contributions. Namely, we perform the perturbative expansion of the n-point functions of g ζ(t f , x a ) and δ g γ ij (t f , x a ) with a = 1, · · · , n, using the canonical variablesΦ. In this subsection, we adopt the Schrödinger picture. Since all the operators will be in Schrödinger picture, they do not have time dependence. Introducing the unitary operator of the time evolution the n-point functions are expressed as Here, we introduce the eigen state of ζ and δγ ij , |ζ c , δγ c . For given values of s and S ij , |ζ c , δγ c also becomes the eigen state of gζ and δ gγ . i.e., where ζ c and δγ c denote the eigenvalues of ζ and δγ ij . Here the time dependence of the operators gζ and δ gγ ij appears through W t (x), xgl and S ij , while ζ(x) and δγ ij (x) are time-independent Schrödinger operators. Since xgl and δγ S, ij depend on the value of s(t) and the path for picked up values of S ij (t ′ ) with t ≤ t ′ ≤ t f , the eigenvalues of the operators gζ and δ gγ ij , s (ev) (t) and S where we labeled the discretized time coordinate from the distant past to t f by negative integers and the one from t f to the distant past by positive integers, with t 0 = t f . For the time being, we focus on the n-point functions for a particular time path of ζ c and δγ c , picking up, at each time step, one representative state among the summed up eigen states in the unit operator (4.10) in Eq. (4.11). Namely, we consider the expectation value for all t. Then, using the corresponding values of s and S ij , we introduce the canonical variablesΦ(x) as defined in the preceding section. UsingΦ, we can replace ζ(xgl ) in Eq. (4.12) withζ(x). Notice that in the canonical system withΦ, we should use the unitary operator of time evolution defined by the HamiltonianH(t), which differs from H(t), as Here, s and S ij are different between the forward and backward paths, and hence the new canonical variablesΦ(x) will differ between them. Furthermore, Eqs. (4.8) and (4.9) imply Next, we write down the expression (4.12) in the interaction picture. Using the unitary operator with an appropriate choice of a lower boundary of the t-integration, the Schrödinger picture fieldsζ(x) and δγ ij (x) are related to the interaction picture fieldsζ I (t, x) and δγ ij I (t, x), respectively, as In the interaction picture, we obtain With gζ I (t) and δ gγ ij I (t) defined as gζ Next, we will show that, when we choose the Euclidean vacuum as the initial state, the n-point functions for g ζ(x) and δ g γ ij (x) can be expanded only in terms of the interaction picture fieldsζ I (x) and δγ ij I (x) with the IR suppressing operators R x . While the interaction Hamiltonian densityH I is messy, the IR regularity can be shown just by using the fact that, as is given in Eq. (3.42),H I is expressed only in terms of and also with the parameters,ṡ Notice that the terms in (4.30) are not suppressed by R x and also that the inverse Laplacian ∂ −2 , which arises from N and N i , may decrease the power of k by 1/k 2 , depending on the choice of the boundary conditions.

Interaction picture fields without the IR suppressing operator
We begin with discussing the first term in Eq. If we can simply replace s(t) and S ij (t) with gζ I (t) and δ gγ ij I (t), respectively in the above expression,ζ I (x) − s(t) and δγ ij I (x) − S ij (t) are reduced toζ I (x) − gζ I (t) and δγ ij I (x) − δ gγ ij I (t), which are combinations suppressed by the IR suppressing operator R x . We will show that the distinctive property of the Euclidean vacuum given in Eq. (4.5) allows us to perform this replacement just by adding terms that are composed only of R xζI (x) and R x δγ ij I (x).
To perform the replacement, we notice that the operator |t a ; ζ c , δγ c I I t a ; ζ c , δγ c | is located next to the interaction Hamilto-nianH I (t a ) as where the abbreviation on the right hand side ofH I denotes operators in the past of t a along the closed time path and the one on the left hand side denotes operators in the future of t a . For notational simplicity, we abbreviate the subscript a in the following discussion. Picking up a singleζ I (x) − s(t) or δγ ij I (x) − S ij (t) from the first term ofH I , given in Eq. (4.33), and using Eqs. (4.28) and (4.29), we rewrite each term as where using A(x), we schematically expressed the operators sandwiched betweenζ I (x) − s(t) or δγ ij I (x) − S ij (t) and |t; ζ c , δγ c I . Since A(x) is a part of the Hamiltonian density in Eq. (4.33), it can be expressed solely in terms of the combinations in Eq. (4.30) and the conjugate momentaπ I andπ ij I . Since gζ I (t) commutes withζ I (x) − s(t), δγ ij I (x) − S ij (t), andπ ij I (x), the non-vanishing contributions in [ gζ I (t), A(x)] arise only from the commutator gζ which yields a local function whose Fourier mode is regular in the IR limit. Repeating this procedure, we can rewrite (ζ I (x) − s(t))A(x) solely in terms ofζ The same argument follows for (δγ ij I (x) − S ij (t))A(x). In this way all the interaction picture fields in the first term ofH I are now expressed by R xζI and R x δγ ij I . Next, we consider the second term of the interaction Hamiltonian (3.42) withṡ andṠ ij . When we discretize the time coordinate, the time derivative should be regarded as the difference between the values at two adjacent time steps. We can express the second term in Eq. (3.42) sandwiched by I t a+1 ; ζ c , δγ c | and |t a ; ζ c , δγ c I as (t a , x). Here, we neglect the terms irrelevant in the continuous limit. Using Eqs. (4.28) and (4.29), we can replace s(t a ) and S ij (t a ) with gζ I (t a ) and δ gγ ij I (t a ) placed next to |t a ; ζ c , δγ c I , while s(t a+1 ) and S ij (t a+1 ) with gζ I (t a+1 ) and δ gγ ij I (t a+1 ) next to I t a+1 ; ζ c , δγ c |. For the same reason as in the previous case, the terms coming from the commutator between [π I (x a )x l ∂ mζI (x a ) + · · · ] and gζ I (t a ) or δ gγ ij I (t a ) only give the IR regular contributions. While, the remaining part becomes Similarly, we can replaceṠ ij (t) in the third term of the interaction Hamiltonian (3.42) with δ gγ ij I (t). Notice that gζ which is manifestly expressed in the IR suppressed form, R xζI (x). To make the IR regularity manifest, in the last equality we added 0 = gζ . In a similar manner we can show δ gγ ij I (t) is also in the IR suppressed form.
In this way, we can show that allζ I s and δγ ij I s in the interaction vertices are multiplied by the IR suppressing operator R x . The argument so far proceeds irrespective of the choice of the initial quantum states. Now, we focus on the distinctive property of the Euclidean vacuum given in Eq. (4.5), which states that the initial states chosen by the boundary condition of the Euclidean vacuum are fixed uniquely and independent of which canonical variables are used for quantization. Therefore, requesting the Euclidean vacuum uniquely determines the initial state irrespective of the picked-up particular path of ζ c and δγ c . Therefore, after the above-mentioned replacements, the possible dependence of the n-point functions on the picked-up path remains only in |t; ζ c , δγ c I I t; ζ c , δγ c |, and hence we can remove the decomposition of unity.

Restricting the interaction vertices to the local region
Next, we will address the inverse Laplacian ∂ −2 . If we choose the boundary conditions for ∂ −2 in N and N i appropriately, N and N i with their argument (t, x) in the region O t can be specified by the fluctuations only within O t . In the general solutions of N and N i given in Eqs. (2.36) and (2.37), the residual gauge degrees of freedom are expressed by arbitrary homogeneous solutions of the Laplace equation, G n (x) and (δ i j − ∂ i ∂ −2 ∂ j )G j n (x). We determine the homogeneous solution G n (x) such that the solution in the observable region O t is given by the convolution between the Green function and the source restricted to the local region, i.e., Similarly, using the transverse part of G i n (x), we can determine the boundary conditions for the remaining ∂ −2 so as to shut off the influence from the region far outside of O t . (For a detailed explanation, see Refs. [66,67].) Then, since all the interaction vertices are confined to the neighborhood of O, the operation of the non-local operator ∂ −2 no longer reduces the power law index with respect to k. Thus, when we choose the Euclidean vacuum as the initial states, we can expand the n-point functions for R x g ζ(t f , x) and R x δ g γ ij (t f , x) only in terms of the interaction picture fields R xζI (x) and R x δγ ij I (x). Since R x g ζ(x) and R x δ g γ ij (x) are not invariant under all the residual gauge transformations, their n-point functions can depend on the boundary conditions of N and N i . However, if we calculate n-point functions for the genuinely gauge invariant operator g R, changing the boundary conditions should not affect the result.

V. REGULARITY OF LOOPS
In this section, we will show that when we choose the Euclidean vacuum as the initial state, the n-point functions of R x g ζ and R x δ g γ ij no longer suffer from the IRdiv, IRsec, and SG. The discussion in this section goes almost in parallel with the one in Sec. IV of Ref. [66], where the regularity of the scalar loops is shown. Here, we briefly highlight the discussion, deferring a more detailed discussion to Ref. [66].

A. Euclidean vacuum from the iǫ prescription
In the preceding section, we introduced the Euclidean vacuum, which satisfies the boundary conditions (4.1). Here, following Ref. [66], we show that these conditions lead to the iǫ prescription in the ordinary perturbative description. For our current purpose, the explicit form of the interaction Hamiltonian densityH I is not necessary. We simply note that all the interaction vertices inH I can be formally expressed as where N s and N t are non-negative integers with N s + N t ≥ 3. Here, λ(t) is a dimensionless time-dependent function which can be expressed only in terms of the horizon flow functions. To discriminate different IR suppressing operators R x s, we added a subscript (m s ) or (m t ) to R x . The spatial indices i m t and j m t will be contracted with other indices i m t′ and j m t′ or with indices in R x , which are abbreviated for notational simplicity. In the following, we use the formal expression (5.1) as the interaction vertices.
Since the boundary conditions for the Euclidean vacuum should also hold at the tree level, the asymptotic form of the positive frequency mode function v α k (t) with α = s or t, in the limit η → −∞, should be proportional to e −ikη(t) . Factoring out this time-dependence, we express v α where we introduced as approximate amplitudes of the curvature perturbation and the graviton perturbation. The function f α k (t) satisfies the regular second order differential equation with the boundary condition Since both the differential equation and the boundary condition of f α k (t) are analytic in k for any t, the resulting function should be analytic as well. In fact, f α k (t) does not have any singularity such as a pole on the complex k-plane as a consequence of the boundary conditions of the Euclidean vacuum.
On the other hand, in the limit −kη(t k ) ≪ 1 the function f α k (t) is proportional to A α (t k )/A α (t), where t k is the Hubble crossing time defined by −kη(t k ) = 1, because the curvature and graviton perturbations should be constant in this limit. Hence, the expansion for small k is in general given by By using Eq. (5.2), the Wightman function for the curvature perturbation is given by and the Wightman function for the graviton is given by where in the second equality we assumed that the quantum state is isotropic.
Using the in-in formalism, the n-point functions can be expanded by the Wightman functions G + s (x, x ′ ), G +t ijkl (x, x ′ ), and their complex conjugates. When we impose the boundary conditions of the Euclidean vacuum, we need to start the vertex integrals at η = −∞. Although the integrands of the vertex integrals are infinitely oscillating in the limit η → −∞, the time integration can be made convergent by adding a small imaginary part to the time coordinate, which is nothing but the ordinary iǫ prescription. To see the convergence of the time integration more explicitly, using the formal expression for the interaction vertex (5.1), we first consider the integral for the vertex which is closest to the past infinity η → −∞(1 − iǫ). The interaction picture fieldζ I (x) included in this vertex is contracted withζ I (x m s ) in vertices labelled by m s = 1, 2, · · · , N s , and give the Wightman function G + s (x m s , x). Similarly, the interaction picture field δγ I i m t j m t (x) included in the vertex is contracted with δγ k m t l m t I (x m t ) in vertices labelled by m t = 1, 2, · · · , N t , and give the Wightman function G + t k m t l m t i m t j m t (x m t , x). Then, the vertex integration with N sζ I s and N t δγ ij I s gives where x m α denotes either x m s or x m t . The Euclidean vacuum condition requires the convergence of this integral in the limit η(t i ) → −∞. Since the Wightman functions contain the exponential factor the integral can be made convergent by adding +iǫ to η(t), which is again exactly what is known as the iǫ prescription. Here, k m s denotes the momentum of G + s (x m s , x) and k m t denotes the momentum of G + t k m t l m t i m t j m t (x m t , x). The vertex integration next to the closest to the past infinity can be done in a similar manner, where N s′ and N t ′ are numbers of scalar and graviton propagators that connect between this second vertex and the vertices other than the first one. If we perform the integration over the time coordinate of the first vertex t up to t ′ , the exponential factor in G + s (x m s , x) or G + t k m t l m t i m t j m t (x m t , x) can be replaced as e ik m α (η(t)−η(tm)) → e ik m α (η(t ′ )−η(tm)) .

(5.11)
Therefore all the Wightman functions connecting the vertices at t ′ or in the past of t ′ with the vertices in the future of t ′ give an exponential factor which is suppressed by adding +iǫ to η(t ′ ). (Here, we mean the future and past in the chronological sense, and not those in the sense of the CTP.) This is again consistent with the boundary condition of the Euclidean vacuum. The same argument can be made for the other vertices as well.
In this subsection, considering the time integration at vertices with fixed momenta of the Wightman propagators, we showed that the boundary condition of the Euclidean vacuum can be imposed in perturbative expansion by employing the iǫ prescription. However, as we will describe in the next subsection, in our proof of the IR regularity, we will perform the momentum integration of the propagator ahead of the vertex integration.

B. IR/UV suppressed Wightman function
Since allζ I (x)s and δγ ij I (x)s in the interaction Hamiltonian are multiplied by the IR suppressing operators R x , the n-point function of R x g ζ(x) and R x δ g γ ij (x) can be expanded by the Wightman function R x ′ ) and their complex conjugates. In this subsection, we calculate the Wightman functions multiplied by the IR suppressing operator, After integration over the angular part of the momentum, the Wightman where we introduced The Wightman function R x R x ′ G + t ijkl (x, x ′ ) is given by a similar expression as Here, before we integrate over the angular coordinates, we replaced the projection tensor P ij with the derivative form: which commutes with R x . We first show the regularity of the k integration in Eqs. (5.12) and (5.13). The regularity of the Wightman function G + s (x, x ′ ) is shown in Ref. [66]. We will see that the same argument also leads to the regularity of G + t ijkl (x, x ′ ).
Since the functions f α k (t) with α = s, t are not singular, the regularity can be verified if the integration converges both in the IR and UV limits. The regularity in the IR limit is guaranteed by the presence of the IR suppressing operator. The IR suppressing operators R x add at least one extra factor of k|η(t)| or eliminate the leading t-independent term in the IR limit, and yield and where we have introduced the spectral indices n s and n t as Thus, the operation of R x makes the k integration in Eqs. (5.12) and (5.13) regular in the IR limit, ensuring the IR regularity. Next, we consider the convergence in the UV limit. When we choose the Euclidean vacuum, the iǫ prescription facilitates the regularization of the UV modes in Eqs. (5.12) and (5.13), because adding a small imaginary part to all the time coordinates as η → η × (1 − iǫ) leads to the replacement Then, the manifest exponential suppression factor is introduced for large k. This UV regulator makes the integral finite for the large k contribution except for the case σ ± (x, x ′ ) = 0, where x and x ′ are mutually light-like. Since the expression of the Wightman functions obtained after the k integration is independent of the value of ǫ, this regulator makes the UV contributions convergent even after ǫ is sent to zero. For σ ± (x, x ′ ) = 0, the integral becomes divergent in the limit ǫ → 0, but the divergence related to the behavior of the Wightman functions in this limit is to be interpreted as the ordinary UV divergences, whose contribution to the vertex integrals must be renormalized by introducing local counter terms. Thus, the are now shown to be regular functions. Since the amplitudes of the Wightman functions with the IR suppressing operator are bounded from above, we can show the regularity of the n-point functions, if the non-vanishing support of the integrands of the vertex integrals is effectively restricted to a finite spacetime region. Since the causality has been established with the aid of the residual gauge degrees of freedom (see Sec. IV B 2), the question to address is whether contributions from the distant past are shut off or not. In short, this question can be rephrased as whether the SG due to the time integral exists or not. To address such a long-term correlation, we discuss the asymptotic behavior of the Wightman functions R , sending t ′ to a distant past. Recall that when σ ± (x, x ′ ) = 0, we can rotate the integration contour with respect to k even toward the direction parallel to the imaginary axis, making ǫ finite. Rotating the direction of the path appropriately depending on the sign of σ ± (x, x ′ ), the integrand shows an exponential decay for k , except for the region where the two points are mutually light-like (see Ref. [66] regarding the estimation of the contribution from this region). The rotation of the k integration contour can be done without hitting any singularity in the complex k-plane, because the functions f α k (t) are guaranteed to be analytic by construction. If we choose other vacua, this operation yields extra contributions from singularities. After the rotation, the integrations of k on the right-hand sides of Eqs. (5.12) and (5.13) are totally dominated by the wavenumbers with k 1/|η(t ′ )| ≪ 1/|η(t)|. Using Eq. (5.15) which gives the asymptotic expansion in the limit k|η(t)| ≪ 1, we obtain 19) where in the second equality we performed the k integration, rotating the integration contour. Similarly, using Eq. (5.16), we obtain We should emphasize that we did not employ the long wavelength approximation regarding the Hubble scale at t ′ to properly evaluate the modes k of O(1/|η(t ′ )|) as well.

C. Secular growth (SG) due to the time integral
In this subsection, focusing on the long-term correlation, we discuss the convergence of the vertex integrals of the n-point functions for the Euclidean vacuum. We start with the integration of the n-point interaction vertex which is the closest to η = −∞(1 − iǫ). By inserting the expression of the Wightman functions R t ≫ t ′ , given in Eqs. (5.19) and (5.20) into Eq. (5.9), the vertex integral V (1) can be estimated as As we have explained in Sec. IV B 2, the interaction vertices are confined within the observable region, i.e., the non-vanishing support of the integrand is bounded by |x| L t , where L t can be approximated by |η(t)| in the distant past. Thus, we obtain Since we have performed momentum integral first, the exponential suppression for large |η|, required for the Euclidean vacuum, no longer remains. However, picking up η-dependence of the integrand of Eq. (5.22), we still find that the contribution from the distant past is suppressed if When this condition is satisfied, the time integral converges, and the amplitude of V Then, when one of the Wightman propagators is connected to a vertex located in the future of x ′ , i.e., t m > t ′ , the t-integration yields the suppression factor {η(t m s )/η(t ′ )} . We denote the numbers of such scalar and graviton propagators byÑ s andÑ t , respectively.
Similarly, we can evaluate the amplitude of V (2) as Extracting the η ′ -dependent part in the above expression, we obtain Now the generalization proceeds in a straightforward way. For the N v -th vertex, the temporal integration becomes where N s f and N t f , respectively, denote the numbers ofζ I s and δγ ij I s contained in the vertices which have been integrated before the N v -th vertex, M s and M t denote the numbers of the Wightman propagators connected to a vertex with η > η Nv , and λ is the product of all the interaction coefficients contained in the integrated vertices. Thus, the convergence condition is derived as As a simple example, we consider the case where ε 1 is constant. In this case,λ is expressed only in terms of ε 1 and takes a constant value. By assuming M = 1 and using n s − 1 = −2ε 1 , the convergence condition yields 29) where N ≡ N s f + N t f − 2N v and M ≡ M s + M t . In the slow roll limit ε 1 ≪ 1, the above condition is recast into Since all interaction vertices contain at least one propagator which is connected to a vertex in their future, M should be M ≥ 1. Therefore, unless an extremely higher order in perturbation with N > O(1/ε 1 ) is concerned, the contributions from the distant past are suppressed and hence the time integrals at the interaction vertices do not yield the SG. The presence of the above suppression can be intuitively understood in the same way as in the discussion for the loops of the curvature perturbation [66]. When we choose the Euclidean vacuum as the initial state, both the IR and UV modes in the Wightman functions are suppressed and then only the contributions around the Hubble scale at each time are left unsuppressed. Being affected only by the modes around the Hubble scale, i.e., k|η| ≃ k/e ρρ = O(1), the Wightman functions R ijkl (x, x ′ ) are necessarily suppressed when η(t)/η(t ′ ) ≪ 1. This is because if the spacetime points x and x ′ are largely separated in time, any Fourier mode in the Wightman function cannot be of order of the Hubble scale simultaneously at t and t ′ . When we consider the contribution of vertices located far in the past, at least one Wightman function should satisfy η(t)/η(t ′ ) ≪ 1, and therefore it is suppressed. Equation (5.28) shows such suppression by M s scalar propagators and M t graviton propagators. As is shown in Eq. (5.28), as we increase the number of operators included in or connected to the interaction vertex, denoted by N s f and N t f , the contributions from the distant past come to less suppressed. On the other hand, as we increase the number of propagators connected to the vertices around the observation time, labelled by M s and M t , the contributions from the distant past are more suppressed. When N is sufficiently large, i.e., N > O(M/ε 1 ), the suppression due to M propagators can be overwhelmed by the large amplitude of the fluctuation, which increases when the energy scale of inflation increases as in the far past. However, we should also stress that the SG never appears in the slow roll inflation, unless the order of perturbative expansion N takes an extremely large value such as 1/ε 1 ≃ O(10 2 ).

VI. CONCLUDING REMARKS
In this paper, we addressed the regularity of the graviton loops. We showed that when we choose the Euclidean vacuum as the initial state, similarly to the curvature perturbation, the graviton perturbation does not cause the IRdiv and IRsec in the n-point functions of genuine gauge invariant operators. In the absence of the graviton, simply performing the dilatation transformation provides the new set of canonical variables in which all ζs in the Hamiltonian are shifted by the free parameter s. Presence of this new set of canonical variables is important to show the IR regularity for the Euclidean vacuum. Extending this previous result to the graviton perturbation, we provided the new set of canonical variables whose Hamiltonian includes the curvature perturbation and the graviton perturbation with the shifts by arbitrary time-dependent parameters s and S ij , respectively. Then, following a similar argument to the one in Ref. [66], we established the IR regularity, i.e., the absence of the IRdiv and IRsec for the Euclidean vacuum to any order of perturbation. We also showed the absence of the SG in the slow roll inflation, at least, unless the extremely high orders in perturbation are concerned.
As is argued also in Ref. [66], when we evaluate the SG, considering only the superH modes is not sufficient, because all modes are subH modes when we send the initial time t i to the past infinity. In Sec. V C, to evaluate the SG, we used the Wightman functions obtained in Sec. V B. These Wightman functions R ijkl (x, x ′ ) are shown to take finite values, as long as the two arguments x and x ′ are not mutually light-like.In this paper and also in Ref. [66], assuming that these UV divergences will be renormalized by local counter terms, we did not explicitly examine the contributions from the singular UV modes. We leave a detailed discussion about the UV renormalization for future study. (See Refs. [76,108], where the UV regularization is discussed. ) In this paper, we considered the inflationary universe as the background spacetime. When we take an exact de Sitter space as the background spacetime without introducing a scalar field, the curvature perturbation will disappear, while the graviton perturbation can still exist. Even if we take the de Sitter limit and pure gravity case is concerned, the accumulation of residual gauge degrees of freedom is still an issue of debate. It has been claimed that the IR graviton can become a trigger of the running of the coupling constant. For instance, in Ref. [10], Tsamis and Woodard claimed that the IR graviton can screen the cosmological constant, suggesting a possibility that the cosmological constant problem might be dynamically solved. In our forthcoming publication, we will address the IR issues of the graviton in the exact de Sitter background and discuss whether the screening of the cosmological constant can still exist even if we request the genuine gauge invariance.
Finally, we make several comments on the quantum states allowed from the IR regularity conditions. We have seen that when we choose the Euclidean vacuum as the initial state, the n-point functions of the genuine gauge invariant operator become IR regular. Then, the question is whether the regularity can be maintained for other initial states or not. In the simple setup adopted in Appendix A, which immediately ensures the standard commutation relations for the interaction picture fields, we found that requesting the IR regularity of the graviton loops yields the same condition on the mode function v s k that was requested from the IR regularity of the loop corrections due to the curvature perturbation. (In Ref. [63], we claimed that the IR regularity of the graviton loops does not yield any condition on v s k . However, as is mentioned in Appendix A, in Ref. [63], we chose an alternative heuristic iteration scheme which does not immediately guarantee the standard commutation relations for the interaction picture fields. Therefore there is no contradiction with the current result.) It will be intriguing to elaborate how strictly the IR regularity condition constrains the quantum state in the inflationary universe. We will also leave this issue for future study. (See also the studies for the scalar field by Einhorn and Larsen in Refs. [109,110] and by Marolf et al. in Ref. [111].) Here, we compute the two-point function of R x g ζ(x) by solving the Heisenberg operator equations of motion for ζ and δγ ij . Using the retarded Green functions G R (x, x ′ ) and G R ijkl (x, x ′ ) given by we can solve the equations of motion for ζ and δγ ij , employing the initial conditions (A1) and (A2) as with where the explicit form of the non-linear source terms S NL (x) and S NL ij (x ′ ) will be derived later. Evaluating Eqs. (A9) and (A10) iteratively, we can obtain expressions for ζ and δγ ij , respectively. Inserting thus obtained solution ζ and δγ ij into Eq. (3.45), we can perturbatively compute g ζ(x) as g ζ(x) = ζ I (x) + g ζ 2 (x) + g ζ 3 (x) + · · · , where g ζ n (x) represents the term that consists of n interaction picture fields. Expanding the interaction picture fields ζ I and δγ ij I as in Eqs. (2.9) and (2.19), the initial vacuum state is defined by Notice that the n-point functions computed by using the Heisenberg operator solved with the retarded Green function can be formally shown to agree with those calculated in the in-in formalism (see, for instance, Appendix of Ref. [65]). Using Eq. (A11), the one-loop contributions to the two-point function of R x g ζ(x) are given by R x1 g ζ(x 1 )R x2 g ζ(x 2 ) 1loop = R x1 g ζ 2 (x 1 )R x2 g ζ 2 (x 2 ) + R x1 ζ I (x 1 )R x2 g ζ 3 (x 2 ) + R x1 g ζ 3 (x 1 )R x2 ζ I (x 2 ) . (A13) to denote the approximate equality neglecting the terms which yield neither ζ2 I nor δγ 2 I at the one-loop level [62,63]. Now, we derive approximate equations of motion for ζ and δγ ij . In the following, we will use with α = s, t. Here, Q I and Q ′ I are either ζ I or δγ ij I and R x is a derivative operator which suppresses the IR modes. Equation (A17) can be proved as follows. The Fourier transformation of L −1 R,s Q I R x Q ′ I is proportional to where Q I k and (RQ ′ I ) k denote the Fourier modes of Q I and RQ ′ I . Since (RQ I ) k−p (t ′ )− (RQ ′ I ) k (t ′ ) is suppressed and Q I,p becomes time-independent in the limit p → 0, the IR relevant piece of the integrand of the momentum integral can be recast into Similarly, we can prove Eq. (A17) also for L −1 R,t . In the following discussion, we will also use the approximate identity In the one-loop corrections to R x g ζ, δγ ij 2 contributes only through g ζ 3 , and δγ ij n with n ≥ 3 do not contribute. Since at least one of two interaction picture fields included in δγ ij 2 is suppressed by R x , we find δγ ij 2 IR ≈ 0. Then, we find and hence the one-loop corrections can be given without computing the non-linear contributions in δγ ij .
Next, we derive an approximate equation of motion for ζ. Under the equality IR ≈, the non-linear action is reduced to S IR ≈ M 2 pl dtd 3 x e 3(ρ+ζ) ε 1 (∂ t ζ) 2 − e −2(ρ+ζ) e −δγ ij ∂ i ζ∂ j ζ , where the terms with more than two fields with differentiation, which give neither ζ2 I nor δγ ij I δγ kl I , are abbreviated. The variation of the above action gives the equation of motion as with L s ≡ ∂ 2 t + (3 + ε 2 )ρ ∂ t − e −2ρ ∂ 2 , and where the last term is added so that the solution satisfies the second condition in Eqs. (A1) [65].
To rewrite the terms with x i ∂ j L −1 R,s in g ζ 3 into a more tractable form, we use the identity which obviously holds if L −1 R,s L s can be replaced with unity. In general, for we have L s δ R = 0, and hence δ R is a homogeneous solution of the second order differential equation i.e., L s δ R = 0. Since δ R and ∂ t δ R are both zero at the initial time, which is automatically satisfied by the definition of the retarded integral L −1 R,s , we can confirm that δ R vanishes for all t ≥ t i . Using L s , x i ∂ j = −2∇ i ∇ j , the right hand side of Eq. (A37) is further rewritten as