%% ****** Start of file apsguide4-1.tex ****** %
%%
%%   This file is part of the APS files in the REVTeX 4.1 distribution.
%%   Version 4.1p of REVTeX, March 2010.
%%
%%   Copyright (c) 2009, 2010 The American Physical Society.
%%
%%   See the REVTeX 4.1 README file for restrictions and more information.
%%
%\documentclass[preprint,secnumarabic,amssymb, nobibnotes, aps, prl]{revtex4-1}
\documentclass[reprint,secnumarabic,amssymb, nobibnotes, aps, prl]{revtex4-1}
%\usepackage{acrofont}%NOTE: Comment out this line for the release version!


%
% \usepackage[usenames,dvipsnames]{pstricks}
% \usepackage{epsfig}
% \usepackage{pst-grad} % For gradients
% \usepackage{pst-plot} % For axes
% % Change this value to rescale the drawing.



\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math
\usepackage{amssymb}
\usepackage{amsmath,amsthm,amsfonts} 
\newcommand{\revtex}{REV\TeX\ }
\newcommand{\classoption}[1]{\texttt{#1}}
\newcommand{\macro}[1]{\texttt{\textbackslash#1}}
\newcommand{\m}[1]{\macro{#1}}
\newcommand{\env}[1]{\texttt{#1}}
\setlength{\textheight}{9.5in}
\def\R{\mathbb{R}}
\def\E{\mathcal{E}}
\renewcommand{\H}{\mathcal{H}}
\def\eqdef{\stackrel{\mbox{\tiny def}}{=}}     % egal a ... par definition
\newcommand{\ket}[1]{|\kern.3ex#1\kern.3ex\rangle}
\newcommand{\bra}[1]{\langle\kern.3ex #1 \kern.3ex|}
\newcommand{\scalar}[2]{\langle\kern.3ex #1 \kern.3ex|\kern.3ex#2\kern.3ex\rangle}
\DeclareMathOperator*{\under}{\longrightarrow} 




\begin{document}

\title{Infrared behavior and gauge artifacts in de Sitter spacetime \\
I. The photon field}

\author{A. Youssef}
\email{ahmed.youssef@apc.univ-paris7.fr}
\affiliation{Universit\'e Paris Diderot Paris 7, Laboratoire APC, B\^atiment Condorcet, Case 7020, 75205 Paris Cedex 13, France}
\date{\today}
\begin{abstract}
We study the infrared (long distance) behavior of the free photon field in de Sitter spacetime. Using a two-parameter family of gauge fixing terms, we show that the infrared (IR) behavior of the two-point function is highly gauge-dependent and ranges from vanishing to growing. This situation is in disagreement with its counterpart  in flat spacetime, where the two-point function vanishes in the IR for any choice of the gauge fixing parameters. A criterion to isolate the ``physical'' part of the two-point function is given and  is shown to lead to a well-behaved two-point function in the IR. 
%We show that in de Sitter space, unlike flat space, the existence of a natural infrared scale, the Hubble scale, leads to an  unusual IR behavior for the massless photon field, namely non-vanishing correlations between largely separated points.  We show that this nothing but a gauge artifact.  
\end{abstract}
%\pacs{04.62.+v}
\maketitle
%\tableofcontents
\section{Introduction}
The study of quantum field theory in the de Sitter background is of paramount importance to the understanding of  the early universe as well as its present accelerated expansion. One of the most striking, yet poorly understood,  aspects of de Sitter QFT is the so called \emph{infrared problem}. The simplest and most serious manifestation of this problem being the non-vanishing of the correlation functions for largely - spacelike and timelike -separated points. This is for instance a well-known fact for the massless minimally coupled (mmc) scalar and the graviton fields. In this letter we study the IR behavior of the photon field in de Sitter space and show that it exhibits similar IR pathologies. More importantly we are able to show that these are purely gauge artifacts. 

The organization of the letter is as follows. First we compute the correlation function using the $R_\lambda$ gauges $\lambda (\nabla_a A^a)^2/2$. Then we generalize the computation to the larger family of higher derivative gauge fixing terms $ \lambda (\nabla_a A^a)^2/2 +\xi R^2 (\nabla_a A^a)\Box (\nabla_a A^a)/2$. Finally we propose a decomposition of the two-point function into a physical and a gauge part and we show that the physical part is vanishing in the IR. 

The $d$-dimensional de Sitter spacetime $X_d$ of ``radius'' $R$ can be identified with the real one-sheeted hyperboloid in the $d+1$ Minkowski spacetime $M_{d+1}$
$$
X_d=\left\lbrace x \in \R^{d+1}, \eta_{\mu \nu}x^\mu x^\nu=R^2\right\rbrace
$$
where $\eta_{\mu \nu}$ is  the ``mostly plus'' flat metric  $\textrm{diag}(-1,+1,\cdots,+1)$. Let $\mu(x,x')$ denote the geodesic distance between the points $x$ and $x'$ and $g_{ab}$ be the d-dimensional de Sitter (or d-dimensional Minkowski if $R\to \infty$) metric. 
We will compute the Wightman two-point function $W_{aa'}(x,x')=\langle\Omega\vert A(x)A(x')\vert \Omega \rangle$ in different gauges. In the previous formula $\vert \Omega \rangle$ designates the Euclidean (Bunch-Davies) vacuum state. This vacuum being de Sitter invariant, $W_{aa'}(x,x')$ is a maximally symmetric bi-tensor. As shown in detail in \cite{Allen86}, all such bi-tensors can be expressed in terms of the parallel propagator $g_{aa'}(x,x')$ between $x$ and  $x'$, and the two vectors $n_a=\nabla_a \mu$ and $n_{a'}=\nabla_{a'} \mu$. Our two-point function is thus written as
$$
W_{aa'}(x,x')=\alpha_1(\mu) g_{aa'}+\alpha_2(\mu) n_a n_{a'} 
$$
where, for spacelike separations ($\mu^2>0$),  $\alpha_1$ and $\alpha_2$ are scalar functions of the geodesic distance $\mu$ only. We will frequently use the following formulas \cite{Allen86}:
\begin{align*}
g_a^{\:\:b'} n_{b'}&=-n_{a'},  \qquad  \qquad   \qquad \:g_{ab}=g_a^{\:\:c'}g_{c' b}\\
\nabla_a n_b&=A (g_{ab}-n_a n_b), \quad \nabla_a n_{b'}=C (g_{ab'}+n_a n_{b'}) \\
 \nabla_a g_{bc'}&=-(A+C) (g_{ab}n_{c'}+g_{ac'}n_b)
 \end{align*}
where 
$$
 A=\frac{\cot(\mu/R)}{R},  \quad \text{and} \quad C=-\frac{\csc(\mu/R)}{R}.
$$
%$$
%\begin{array}{cccc}
%g_a^{\:\:b'} n_{b'}&=-n_{a'},  &g_{ab}&=g_a^{\:\:c'}g_{c' b}\\
%\nabla_a n_b&=A (g_{ab}-n_a n_b), &\nabla_a n_{b'}&=C (g_{ab'}+n_a n_{b'}) \\
% &\nabla_a g_{bc'}&=-(A+C) (g_{ab}n_{c'}+g_{ac'}n_b)&\\
% A&=\frac{\cot(\mu/R)}{R}, &C&=-\frac{\csc(\mu/R)}{R}
%\end{array}
%$$
It is also useful to introduce the quantity $z$ given for spacelike separations ($\mu^2>0$ or $0<z<1$) by
\begin{equation}
\label{defz}
z=\cos^2 \left( \frac{\mu}{2R}\right).
\end{equation}
All of our calculations are performed in the spacelike region $0<z<1$. They can be extended further by analytic  continuation in the $z$ variable to the cut complex plane $\mathbb{C} \setminus (1,\infty)$. Finally, the Feynman propagator $G^F_{a a'}=i \langle \Omega \vert T A_a(x) A_{a'}(x')\vert \Omega \rangle $, can be obtained as the limiting value of $W_{aa'}$ when approaching the branch cut $z>1$ from above.

We end our introduction by noting that several, quite interesting, recent works (see \cite{MarolfHollands} among many others) were devoted to the study of infrared effects in de Sitter. In these articles however, only massive (although \emph{interacting}) fields are considered. IR pathologies are of course much stronger for massless fields and they appear already at the linear level.
\section{The $R_\lambda$ gauges}
The gauge fixed action describing the massless vector field $A_a$ is
\begin{equation}
\label{action1}
S_\lambda=\int d^4x \sqrt{-g} \left( \frac{1}{4} F_{ab}F^{ab}+ \frac{\lambda}{2} (\nabla_a A^a)^2\right)
\end{equation}
where $F_{ab}=\nabla_a A_b-\nabla_b A_a$. The resulting equation of motion is $D_{a b} A^b=0$ where 
\begin{equation}
\label{operatorD}
D_{a b}= g_{ab} \Box-\frac{(d-1)}{R^2} g_{ab}+(\lambda-1) \nabla_a \nabla_b.
\end{equation}
\paragraph{\textbf{Flat space.}}
In flat space ($R \to \infty$), the Fourier transform of the Feynman propagator reads 
$$
G^F_{aa'}(k)=\frac{1}{k^2-i\epsilon} g_{aa'}+\frac{1-\lambda}{\lambda} \frac{1}{(k^2-i\epsilon)^2}k_a k_{a'}.
$$
Using the (massless limit) of the scalar Feynman propagator (given, as always, for spacelike separations) 
\begin{equation}
\label{scalargreen}
G^F(x,x')=\int \frac{d^4k}{(2\pi)^4} \frac{e^{i k (x-x')}}{k^2+m^2-i \epsilon}=\frac{i m}{(2\pi)^2} \frac{K_1(m \mu)}{\mu}
\end{equation}
(where $K$ is the modified Bessel function of the second kind) and its derivatives with respect to $x$ and $x'$ we get the coordinate-space two-point function
\begin{equation}
W_{aa'}= \left[\frac{\lambda +1}{8 \pi ^2 \lambda  \mu ^2}\right] g_{aa'}-\left[\frac{\lambda -1}{4 \pi ^2 \lambda  \mu ^2}\right] n_{a} n_{a'}.
\end{equation}
For all values of $\lambda$, this expression vanishes for large timelike and spacelike separations as expected.

\paragraph{\textbf{De Sitter space.}}
Following \cite{Allen86}, we write the two-point function in de Sitter space as 
\begin{equation}
\label{generalW}
W^{[\lambda]}_{a a'}(x,x')= \alpha_1^{[\lambda]}(z) \: g_{a a'}+\alpha_2^{[\lambda]}(z) \: n_a n_{a'}. 
\end{equation}
where $\alpha_1,\alpha_2$ are scalar functions of the invariant quantity $z$ defined in \eqref{defz}. The equation of motion $D_{ab} W^{b}_{\:\:a'}=0$ implies two independent equations on $\alpha_1$ and $\alpha_2$. These equations are most easily solved in the Feynman gauge $\lambda=1$ and the result is found in \cite{Allen86} to be
\begin{align*}
\alpha_1^{[\lambda=1]}(z)&=\frac{1}{48 \pi^2 R^2} \left[ \frac{3}{1-z} +\frac{1}{z}
+\left( \frac{2}{z}+\frac{1}{z^2}\right) \ln(1-z)\right]\\
\alpha_2^{[\lambda=1]}(z)&=\frac{1}{24 \pi^2 R^2 } \left[ 1-\frac{1}{z} +\left( \frac{1}{z}-\frac{1}{z^2} \right)\ln(1-z)\right]
\end{align*}
Instead of vanishing for largely separated points, the correlation goes to a non-zero constant function
\begin{equation}
\label{IR1}
\lim_{z\to \infty} W^{[\lambda=1]}_{aa'}(x,x')=\frac{1}{24 \pi^2 R^2 } \: n_a n_{a'}.
\end{equation}
This constant vanishes in the flat spacetime limit ($R\to \infty$) as expected. The non-vanishing of the correlator  in de Sitter space seems to contradict our experience with the more familiar flat space QFT. It is indeed a general result of the latter, sometimes referred to as the cluster decomposition principle \cite{Strocchi}, that the correlation functions must decay in the IR. To our knowledge, only very special QFT models, e.g confinement models like the massless scalar in two-dimensional flat space, seem to escape this result.

As explained in the introduction, it turns out that this IR pathological behavior of the photon field is typical of de Sitter massless fields. The most notable examples are the mmc scalar and the graviton fields. The correlation functions of these fields exhibit even stronger pathologies as they grow in the IR (see for instance \cite{Moschella,Staro79,Antoniadis91,Higuchi}). 

This ``IR problem'' is at the heart of an open discussion in the community. It led some authors \cite{Antoniadis07} to speak of a  ``striking violation of cluster decomposition properties of the de Sitter invariant vacuum state''. It has also been shown that this growing of the graviton two-point function render scattering amplitudes divergent \cite{ilio87}. As a consequence, the authors of \cite{Antoniadis07}, among many others, suspect  a \emph{quantum instability} of de Sitter space, meaning that the de Sitter geometry is not a stable ground state of quantum gravity with a cosmological constant. 

Other workers in the field, mainly motivated by the fact that the correlations of \emph{some} gauge invariant quantities fall-off in the IR limit, claim that the non-vanishing of the two-point functions in the IR is nothing but a gauge artifact. For instance, the electromagnetic field correlation function is well-behaved in the IR \cite{Allen86}:
$$
\langle F_{ab} F^{a' b'}\rangle = \frac{1}{8 \pi^2 R^4} \frac{1}{(1-z)^2}  \left( g_{[a}^{\:\:[a'}g_{b]}^{\:\:b']} + 4\: n_{[a} g_{b]}^{\:\:[b'} n^{a']} \right).
$$
The same is true for the graviton field as there is no growing terms in the correlation
function of the Riemann tensor $\langle R_{abcd}(x) R_{a'b'c'd'}(x')\rangle$. 

In the following we will study in detail the IR problem in the case of the photon field in de Sitter. The IR behavior of this field received considerably less attention in the literature than its counterpart for the mmc and the graviton field. This is certainly due to the importance of these two-fields in inflation theory. It is also probably due to the fact that in actual QED calculations, the non vanishing term \eqref{IR1} will not contribute to scattering processes (see last section).  This model will however prove to be precious in realizing that - at least for some situations - serious IR pathologies in de Sitter space are - without doubt - gauge artifacts.  

%First we show that it is indeed possible to  cure the IR problem of the photon field at the level of the two-point function. 

The equation of motion \eqref{operatorD} applied to the two-point function \eqref{generalW} leads after some algebra to a coupled system of differential equations. Hereafter we will only need one of these equations, namely:
\begin{align}
\label{eom1}
\nonumber &\mathbf{H}\left[1,d-2,\frac{1}{2}+\frac{1-\lambda}{2}\right] \alpha_1+\frac{\lambda -1}{2}  \alpha_2' \\
&+\frac{\lambda -\lambda _d}{\lambda_d-1} \left[ \frac{\alpha _1 }{z}+ \frac{1}{2}\left( \frac{1}{ 1-z}-\frac{1}{ z}\right)\alpha_2\right]=0 
% \\ \nonumber&\mathbf{H}\left[ 0,d-1,\frac{d}{2}+\frac{\lambda-1}{2}\right] \alpha_2-\left[ -\frac{(\lambda -1) (d+\lambda -1)}{2 \lambda }\right] \alpha_1'\\
%\nonumber&+\left[\frac{(d-1) (2 z-1) \lambda ^2+2 \lambda  \left[(3-d)z+d-2\right]+d-1}{4 \lambda 
%   (z-1) z }\right] \alpha_2\\
%\label{eom2}
%&-\left[ \frac{2 (d-2) (z-1) \lambda +(d-1) \lambda ^2-d+1}{2  \lambda z }\right] \alpha_1=0
\end{align}
where the prime denote derivation with respect to the variable $z$, $\lambda_d=\dfrac{d-3}{d-1}$ and  $\mathbf{H}$ is the hypergeometric operator defined by
$$
\mathbf{H}[a,b,c]= z(1-z)\frac{d^2}{dz^2}+\left[ c-(a+b+1) z\right]\frac{d}{dz}-a b 
$$
%$
%\quad \lambda_d'=2\frac{d-2}{d-1}.
%$
Another, moreover first order, relation between $\alpha_1$ and $\alpha_2$ is obtained by considering the two-point function of the field strength
$$
\langle F_{ab} F^{a' b'}\rangle = \beta_1 \: g_{[a}^{\:\:[a'}g_{b]}^{\:\:b']} + \beta_2 \: n_{[a} g_{b]}^{\:\:[b'} n^{a']}.
$$
In fact, the equation of motion (in vacuum)  $\nabla^a F_{ab}=0$ implies the relation \cite{Allen86}
%\begin{equation}
%\label{relation}
%z(1-z) \left[ \frac{1}{2} \beta_2'-\beta_1' \right]-(d-2) \left[ (1-z) \beta_1+\frac{2z-1}{4} \beta_2\right]=0
%\end{equation}
\begin{align}
\label{relation}
%\mathbf{H}\left[d-1,2,1+\frac{d}{2}\right] \beta_1&=0\\
\alpha_2&=\frac{z}{8 \pi ^2 R^2 (z-1)}+2 (1-z) \left(z \alpha_1'+\alpha _1\right)
\end{align}
We now specialize our results to the $d=4$ case for simplicity. The general solution in the Euclidean vacuum is found to be  
\begin{align*}
\label{solution1}
\alpha_1^{[\lambda]}(z)=\frac{1}{48 \pi^2 R^2 \lambda} &\left[ \frac{3(\lambda+1)}{2(1-z)} +\frac{3\lambda-1}{2z}+\right.\\
 & \left.(3\lambda-1)\left( \frac{1}{z}+\frac{1}{2 z^2}\right) \ln(1-z)\right]\\
\alpha_2^{[\lambda]}(z)=\frac{1}{24 \pi^2 R^2 \lambda} &\left[ 1-\frac{3(\lambda-1)}{2(1-z)} -\frac{3\lambda-1}{2z}+\right.\\
&\left.(3\lambda-1)\left( \frac{1}{2z}-\frac{1}{2 z^2} \right)\ln(1-z)\right]
\end{align*}
and the IR behavior is given by
\begin{equation}
\lim_{z\to \pm \infty} W^{[\lambda]}_{aa'}(x,x')=\frac{1}{24 \pi^2 R^2 \lambda} \: n_a n_{a'}.
\end{equation}
This result means that one can cure any IR bad behavior by going to the \emph{Landau gauge} $\lambda \to \infty$. More importantly, we see that the IR behavior of the correlation function is \emph{pure gauge}, thus asking whether or not the two-point function is well or ill-behaved in the IR is in itself -at least partially-a misleading question.  

Finally we note that there exists a special choice of the gauge fixing parameter, namely $\lambda=\frac{1}{3}$ (more generally $\lambda_d$ in $d$ dimensions), that cancels logarithmic terms and gives a particularly simple two-point function. We note, perhaps as a curiosity for the time being, that for large dimensions $d$, $\lambda_d$ goes to the Feynman gauge $\lambda=1$, which is the simplest gauge in flat space. 
%$$
%W^{[\lambda=\frac{1}{3}]}_{aa'}(x,x')=\frac{1}{8\pi^2 R^2} \left[ \frac{1}{1-z} \left(g_{aa'}+n_a n_{a'}\right)+ n_a n_{a'}\right]
%$$
%It is interesting that in the limit $d\to \infty$, $\lambda_d$ goes to the Feynman value $\lambda=1$. This is indeed expected: as the dimension increase, the interaction gets diluted. 
%Finally we comment on the flat space situation. 
%When considering higher loops diagrams, The use of the ill-behaved two-point function can be catastrophic?.

\section{More pathological gauges}
%In the higher derivative gauge, the two-point function even grows. This is not the case in flat space. 
%Given the technical complexity in solving exactly the systems of differential equations on the $\alpha_i$, specially for more exotic gauge choices.
We now consider a larger family of gauge fixing terms. Our motivation is that the IR pathology  we exhibited in the last section, namely that the correlation function tends to a constant, might seem -I believe mistakenly- mild enough not to worry about. We will show that in this more general gauge, the two-point function is actually growing in the IR, exactly like the mmc scalar or the graviton fields. 
%Second, it will motivate our proposal to use the Landau gauge in a more systematic manner as we will shortly explain.
 We consider the action 
\begin{align}
\label{action2}
S_{\lambda,\xi}=S_\lambda+\frac{\xi}{2 m^2} \int d^4x \sqrt{-g}  (\nabla_a A^a)\Box (\nabla_a A^a)
\end{align}
where we have added a higher derivative gauge fixing term. Similar higher derivative gauge fixing are occasionally used in QCD and the electro-weak theories. They were also used to study flat space perturbative quantum gravity \cite{ilio86}. We note that if we want $\lambda$ and $\xi$ to be dimensionless, such higher derivative gauge fixing terms require the introduction of a mass quantity $m$ . In flat space $m$ is introduced by hand, while in de Sitter space the inverse of the de Sitter radius plays naturally this role:  $m^2=\frac{1}{R^2}$. This observation renders the introduction of higher derivative gauge fixing somehow more natural in the de Sitter case. The equation of motion reads
\begin{equation}
\label{eom2}
\left[\Box g_{ab}+(\lambda-1) \nabla_a \nabla_b +\frac{\xi}{m^2} \Box  \nabla_a \nabla_b \right]A^b=0
\end{equation}
\paragraph{\textbf{Flat space.}}
The flat space Feynman propagator is found in Fourier space to be
 \begin{align*}
\tilde{G}^F_{a a'}(k)=&\frac{ g_{aa'}}{k^2-i\epsilon}+\Bigg[\frac{1-\lambda }{\lambda\:  (k^2-i\epsilon)^2 } \\
&+\frac{\xi}{m^2 \lambda ^2} \left( \frac{1}{k^2-i\epsilon}- \frac{1}{k^2-\frac{m^2 \lambda}{\xi}-i\epsilon}\right)\Bigg] k_a k_{a'}
\end{align*}
The reader will note the relative negative sign  typical of higher derivative theories in front of the last propagator.
Using \eqref{scalargreen} and its derivatives the coordinate-space two-point function is then found to be 
\begin{align*}
W_{aa'}&= \left[ \frac{1}{2 \pi ^2 m^2 \lambda  \mu ^4 \nu ^2}+\frac{1}{2 \pi ^2 m^2 \lambda  \mu ^4 \nu ^2}+\frac{K_2(i m \mu  \nu )}{4 \pi ^2 \lambda  \mu ^2}\right] g_{aa'}+\\
   &\left[\frac{2}{\pi ^2 m^2 \lambda  \mu ^4 \nu ^2}+\frac{1-\lambda }{4 \pi ^2 \lambda  \mu ^2}+\frac{K_0(i m \mu  \nu )}{\pi ^2 \lambda  \mu ^2}\right. \\
   &\left. +i \left(\frac{ m \nu }{4 \pi ^2 \lambda  \mu }-\frac{2 }{\pi ^2 m
   \lambda  \mu ^3 \nu }\right) K_1(i m \mu  \nu )\right] n_{a} n_{a'}
\end{align*}
where we have introduced the quantity $\nu=(\lambda/\xi)^{1/2}$. For all values of $\lambda$ and $\nu$, this expression vanishes for large timelike and spacelike separations as expected. The short distance singularity is given by
$$
W_{aa'}\sim \frac{1}{8\pi^2 \mu^2} \left( g_{aa'}-2 n_a n_{a'}\right).
$$
We emphasize that the short distance singularity of the two theories \eqref{action1} and \eqref{action2} are different. This is due to the non commutation of the limits $\xi \to 0$ and $\mu \to 0$.
%\paragraph{Maximally symmetric bitensors}
%Following the same methods described above we obtain the following equations on  $\alpha_1$ and $\alpha_2$:
%\begin{align*}
%&\alpha_1^{(4)}+\frac{2(d+1)}{\mu} \alpha_1^{(3)}+\left( m^2\tau^2 +\frac{d^2-1}{\mu^2}\right) \alpha_1''   \\
%&+\left(\frac{(d+1) m^2 \tau^2}{\mu} -\frac{d^2-1}{\mu^3}\right) \alpha_1' 
%=-\frac{(\lambda +1) m^2 \tau^2 \Gamma(d/2)}{2\lambda \pi^{d/2} \mu^d}\\
%& \alpha_2 = -\frac{\mu ^{2-d} \Gamma
%   \left(d/2\right)}{2\pi ^{d/2} } -\mu\: \alpha_1' .
%   \end{align*}
%   where we defined $\tau=(\lambda/\xi)^{1/2}$. The general solution of the homogeneous equation for $\alpha_1$ reads
%\begin{align*}
%C_1 \:\mu ^{-d/2} J_{\frac{d}{2}}(m \tau \: \mu )+C_2\: \mu
%   ^{-d/2} Y_{\frac{d}{2}}(m \tau \: \mu  )+C_3\: \mu
%   ^{-d}+C_4
%\end{align*}
%and a particular solution $\tilde{\alpha}_1$ of the inhomogeneous equation can be found by the method of variation of parameters. For simplicity we now specify our results to $d=4$ and find
%$$
%\tilde{\alpha}_1= \frac{(\lambda +1) \xi }{2 \pi ^2 m^2 \lambda ^2 \mu
%   ^4}+\frac{\lambda +1}{8 \pi ^2 \lambda  \mu ^2}.
% $$
%Using the asymptotic behavior of Bessel functions
%$$
%J_\nu(\tau \mu), Y_\nu(\tau \mu) \approx \frac{2^{-\nu }  \tau ^{\nu }}{\Gamma (\nu +1)} \mu ^{\nu }, \quad \mu \to 0.
%$$
%several boundary conditions can be chosen, notably the following two-point function have a $\mu^{-2}$ leading singularity in the UV.
% \begin{align*}
% \alpha_1(\mu)=& \frac{\frac{8 \pi  c_2 \lambda }{\tau ^2}+\lambda +1}{2 \pi ^2
%   \lambda  \mu ^4 \tau ^2}+\frac{\left(\frac{8 \pi  c_2 \lambda
%   }{\tau ^2}+\lambda +1\right) Y_2(\mu  \tau )}{8 \pi  \lambda  \mu
%   ^2}\\
%   &-\frac{c_1 J_2(\mu  \tau )}{\mu ^2 \tau ^2}+\frac{\lambda
%   +1}{8 \pi ^2 \lambda  \mu ^2}\\
% \alpha_2(\mu)=& \frac{2 \left(8 \pi  c_2 \lambda +\lambda  \tau ^2+\tau
%   ^2\right)}{\pi ^2 \lambda  \mu ^4 \tau ^4}+\frac{1-\lambda }{4
%   \pi ^2 \lambda  \mu ^2}\\
%   & -\frac{c_1 J_3(\mu  \tau )}{\mu  \tau }+ \frac{8 \pi  c_2 \lambda +\lambda  \tau ^2+\tau ^2}{8 \pi  \lambda 
%   \mu  \tau } Y_3(\mu  \tau )
%\end{align*}
% Using the asymptotic form of Bessel functions
% \begin{align*}
% J_\nu(\mu) &\approx \sqrt{\frac{2}{\pi \mu}} \sin \left( \frac{\pi}{4}-\frac{\pi \nu}{2}+\mu\right), \quad \mu \to \infty \\
% Y_\nu(\mu) &\approx \sqrt{\frac{2}{\pi \mu}} \sin \left( -\frac{\pi}{4}-\frac{\pi \nu}{2}+\mu\right), \quad \mu \to \infty 
% \end{align*}
% it is clear that the two-point function vanishes in the IR. 

\paragraph{\textbf{De Sitter space.}}
We now consider the action \eqref{action2} in de Sitter space and thus set $1/m^2=R^2$.
Following the same methods described above we obtain the following equations on  $\alpha_1$:
\begin{align*}
\alpha_1^{(4)}&=\frac{4 (\tau -2) (\tau +3)}{(z-1)^2 z^2} \alpha_1+\frac{3 (2 z-1) \left(\tau ^2+\tau -16\right)}{(z-1)^2 z^2} \alpha_1'\\
&-\frac{12-(z-1) z (\tau -8) (\tau +9)}{(z-1)^2 z^2}\alpha_1''-8 \left(\frac{1}{z}+\frac{1}{z-1}\right) \alpha_1'''
%=&\frac{\left(\tau ^2+\tau -2\right) (2 z (3 \lambda -1)+3)+\lambda  (18-5
%   \tau  (\tau +1))}{32 \pi ^2 R^2(z-1)^4 z^2 \lambda }\\
%   \alpha_2=&\frac{z-16 \pi ^2 R^2 (z-1)^2 \left(z \alpha_1'+\alpha_1\right)}{8 \pi ^2 R^2 (z-1)}
\end{align*}
where we introduced the parameter $\tau=\frac{\sqrt{4 \lambda +9 \xi }-\sqrt{\xi }}{2 \sqrt{\xi }}$. Obtaining this equation requires some elaborate algebraic manipulations. The calculations in this section were thus verified with the tensor algebra system xAct on Mathematica \cite{xAct}. 
%The general solution of the homogeneous equation on $\alpha_1$ reads
%\begin{align*}
%\alpha_1(z)=&C_1 \frac{P_{\tau }^2(2 z-1)}{(z-1) z}+C_2 \frac{Q_{\tau }^2(2 z-1)}{(z-1) z}\\
%&+C_3 \frac{2 z^3-3 z^2-4}{(z-1)^2 z^2}+C_4 \frac{1}{(z-1)^2 z^2}
%\end{align*}
%where $P_{\tau}^{2}$ and $Q_{\nu}^{2}$ are the Legendre functions of the first and second kind respectively. 
%A particular solution of the inhomogeneous equation can be found by the method of variations of parameters:
%\begin{align*}
%&\frac{(2 z+1) (3 \lambda -1) \log (1-z)}{96 \pi ^2 R^2 z^2 \lambda }+\frac{\tau  (\tau +1)}{24 \pi ^2 R^2 z (\tau -1) (\tau +2)}\\
%&-\frac{(9 \lambda -1) \tau ^2+(9 \lambda -1) \tau -33 \lambda +11}{144 \pi ^2 R^2 (z-1) \lambda  (\tau -1) (\tau +2)}\\
%&\frac{\lambda  \left(\tau ^2+\tau +6\right)+\tau ^2+\tau -2}{192 \pi
%   ^2 R^2 z^2 \lambda  \left(\tau ^2+\tau -2\right)}\\
% &  +\frac{\lambda  (9 \tau  (\tau +1)-48)-11 \tau  (\tau +1)+40}{576 \pi
%   ^2 R^2 (z-1)^2 \lambda  (\tau -1) (\tau +2)}
%\end{align*}
%Note that we skip some lengthy manipulations on special functions that are required to put the solution $\alpha_1$ in such a simple form. The reader who wants to verify that this is indeed a solution has jut to check it directly. 
This equation can be solved in closed form. Using the asymptotic formulas in table \ref{table} the solution verifying: i)regularity at $z=0$ and ii)flat space short distance singularity is found.  For non integer values of $\tau$ it reads
%\begin{align*}
%\alpha_1(z)&=-\frac{\cot (\pi  \tau ) P_{\tau }^2(2 z-1)}{32 \pi  R^2 (z-1) z
%   \lambda  \left(\tau ^2+\tau -2\right)}\\
%&+\frac{Q_{\tau }^2(2 z-1)}{16 \pi ^2 R^2 (z-1) z \lambda  \left(\tau ^2+\tau
%   -2\right)}\\
%   &+\frac{3 \lambda -1}{96 \pi ^2 R^2 \lambda } \left[ \frac{(2 z+1) \log (1-z)}{z^2}+\frac{1}{z}\right]\\
%   &+\frac{1}{32 \pi ^2 R^2 \lambda  (\tau -1) (\tau +2)}\\
%   & \left[ \frac{1}{(z-1)^2}-\frac{(\lambda +1) \tau \left(\tau+1 \right)-2 \lambda}{z-1}\right]
%     \end{align*}
\begin{align*}
\alpha^{[\lambda,\tau]}_1(z)&=\frac{ 2 Q_{\tau }^2(2 z-1)-\pi \cot (\pi  \tau ) P_{\tau }^2(2 z-1)}{32 \pi^2  R^2 \lambda  \:(z-1) z  \left(\tau ^2+\tau -2\right)} \\
   &+\frac{3 \lambda -1}{96 \pi ^2 R^2 \lambda } \left[ \frac{(2 z+1) \ln (1-z)}{z^2}+\frac{1}{z}\right]\\
   &+\frac{1}{32 \pi ^2 R^2 \lambda  (\tau -1) (\tau +2)}\\
   & \left[ \frac{1}{(z-1)^2}-\frac{(\lambda +1) \tau \left(\tau+1 \right)-2 \lambda}{z-1}\right]
     \end{align*}
where $P$ and $Q$ are Legendre functions of the first and second kind respectively.  As before, $\alpha_2$ is obtained by \eqref{relation}. The asymptotic behavior in the IR is given by ($\tau>1$):
\begin{align*}
\alpha^{[\lambda,\tau]}_1(z)\sim&-\frac{2^{2 \tau -5}\cot (\pi  \tau ) \Gamma \left(\tau
   +\frac{1}{2}\right)}{\pi ^{3/2} R^2 \lambda  (\tau +2) \Gamma
   (\tau )}  z^{\tau -2} \\ 
   \alpha^{[\lambda,\tau]}_2(z)\sim&\frac{2^{4 \tau -1} (\tau -1) \tau  \cot (\pi  \tau )
   \Gamma \left(\tau +\frac{3}{2}\right) \Gamma \left(\tau
   +\frac{5}{2}\right)}{\pi ^2 R^2 \lambda  (2 \tau +3) \left(2 \tau
   ^2+5 \tau +2\right) \Gamma (2 \tau +2)} z^{\tau -1} 
%   &+\frac{1}{24 \pi ^2 R^2
%   \lambda }
\end{align*}
Thus for $\tau>2$  for instance, the two-point function grows like $z^{\tau-2}$ in the IR. This strong IR growing is purely a gauge artifact since, as proven in the beginning of this letter, no IR pathologies arise in the Landau gauge. 
\squeezetable
\begin{table}[htb!]
  \centering 
\begin{tabular}{|c|c|c|c|}
\hline
% after \\ : \hline or \cline{col1-col2} \cline{col3-col4} ...
&&&\\
   & $z \to 0$ & $z \to 1$ & $z \to \infty$    \\
&&& \\
\hline
\hline
&&&\\
 $P_{_\tau}^{(2)}(2z-1)$ & $ -\dfrac{\sin(\pi \tau)}{\pi z}$ & $  \dfrac{\Gamma(\tau+3)}{2 \:\Gamma(\tau-1)}  (z-1)$ & $ \dfrac{2^{2\tau} \Gamma(\tau+1/2) }{\sqrt{\pi}\:\Gamma(\tau-1)} \: z^\tau$\\ 
 &&&\\
 \hline
 &&&\\
 $Q_{_\tau}^{(2)}(2z-1)$ & $- \dfrac{\cos(\pi \tau)}{2 z}$ & $\dfrac{1}{2(z-1)}$ & $ \dfrac{\sqrt{\pi}\: \Gamma(\tau+3) }{2^{2(\tau+1)} \Gamma(\tau+3/2)} \dfrac{1}{z^{\tau+1}}$\\
 &&&\\
 \hline
\end{tabular}
  \caption{Asymptotics of the Legendre functions near singular points.}\label{table}
\end{table}

%Use the relations
%$$
%Q^m_\tau(z)=(-1)^m (1-z^2)^{m/2} \frac{d^m}{dz^m} Q_\tau(z)
%$$
%\scalebox{1}
%{\begin{pspicture}(0,-2.5379)(6.4958,2.5779)
%\definecolor{color879}{rgb}{1.0,0.2,0.2}
%\rput(0.6758,-1.9021){\psaxes[linewidth=0.036,ticksize=0.0358cm]{->}(0,0)(0,0)(4,4)}
%\usefont{T1}{ptm}{m}{n}
%\rput(1.1458,2.3829){$\tau$}
%\usefont{T1}{ptm}{m}{n}
%\rput(5.6958,-1.5971){$\lambda$}
%\psline[linewidth=0.04cm,linecolor=color879](0.5558,0.0579)(4.5558,0.0579)
%\end{pspicture} }


\section{A physical decomposition}
\label{sectiondecomposition}
The previously described situation, namely that the IR behavior of the two-point function in de Sitter space contains important gauge artifacts makes it natural to look for a decomposition of the two-point function into a physical and a non-physical parts. We show now that such a decomposition exists and is given by the following rewriting of the two-point function
\begin{equation}
\label{vecphysdecomp}
W^{[\lambda]}_{a a'}(x,x')= g_{aa'}\beta_1(z)+\nabla_a \nabla_{a'} \beta_2(z).
\end{equation}
In fact, most of the physics is included in the ``quantum'' action $W_J$ defined by
$$
\exp\left[\frac{i}{\hbar} W_J\right]=\int D[A_a]  \exp \left[\frac{i}{\hbar} \left(S+\int d^4x \sqrt{-g} \: A_a J^a\right)\right]
$$  
where $J_a$ is an external  conserved current $\nabla_a J^a=0$. The tree-level expression of the effective action is given by
$$
W_J=\int dV_x dV_{x'}  J^a(x) W_{a a'}(x,x') J^{a'}(x')
$$
where $dV_x=d^4x \sqrt{-g(x)}$ is the invariant volume element. Since the current $J_a$ is conserved, integration by parts ensures that only $\beta_1$ will contribute to $W_J$ and it thus fully deserves to be referred to as the ``physical'' two-point function.
%The conservation of the current $J_a$ and the special  decomposition \eqref{vecphysdecomp} of the two-point function implies 
%\begin{align*}
%W_J&=-\int dV_x dV_{x'} \: J^a(x) g_{aa'} J^{a'}(x') \beta_1(z) \\
%       &=\int dV_x \: J^a(x) J_{a}(x) \beta_1(z) 
%\end{align*}
%$$
%\nabla_a \nabla_{a'} f=\frac{f'}{2R^2} g_{aa'}+\frac{1-z}{R^2} (f'+z f'') n_a n_{a'}
%$$
An explicit form of the functions $\beta_1$ and $\beta_2$ can be found by solving the system
\begin{align*}
\alpha_1=\beta_1+\frac{1}{2R^2} \beta_2', \quad R^2 \alpha_2=(1-z)\beta_2' +z(1-z)\beta_2''.
\end{align*}
The solutions are unique if one requires regularity near $z=0$ and we obtain
%\begin{align*}
%\beta_1^{[\lambda]}(z)&=\beta_1(z)=\frac{1}{1-z}+\frac{\ln(1-z)}{z} \\
%\beta_2^{[\lambda]}(z)&= %\end{align*}
\begin{equation}
\label{phys}
W^{\textrm{Phys}}_{a a'}(x,x')=\frac{1}{16 \pi^2 R^2} \left(\frac{1}{1-z}+\frac{\ln(1-z)}{z}\right) g_{a a'}
\end{equation} 
It is a remarkable (and expected) fact that $W^{\textrm{Phys}}_{a a'}$ is $(\lambda,\xi)$ independent and is well behaved in the IR.
In flat space-time, the Feynman gauge $(\lambda=1,\xi=0)$ gives precisely this physical two-point function. The situation in de Sitter is quite different and no choice of the gauge fixing parameters gives the physical two-point function. A natural question is wether a more general gauge fixing term gives directly $W^{\textrm{Phys}}_{a a'}(x,x')$ in de Sitter space. 
\section{Outlook}
We showed that the IR pathological behavior of the photon's two-point function is nothing but a gauge artifact. We suggest that the situation is similar for the graviton field. The IR behavior of the latter is however necessarily  more intricate since it has been shown that the (transverse-traceless) growing term of the two-point function does not cancel out from scattering amplitudes \cite{ilio87,Antoniadis91}. This implies in particular that a decomposition like \eqref{phys} does not lead to a well-behaved ``physical''  two-point function. The IR behavior of the graviton field  will be examined  in \cite{Ahmed}.   

%Is answer is positive. This gauge is a higher derivative one and reads
%\begin{align*}
%S_{\textrm{gf}}=&\frac{1}{2} \left(\nabla {}_{a}A{}^{a}\right)^2 - 4 A{}_{a} \Box A{}^{a} - \frac{R^2}{2}  \left(\nabla 
%{}_{a}A{}^{b}\right) \Box \left(\nabla {}^{a}A{}_{b}\right) + \\
%&\frac{R^2}{2}  \left(\nabla {}_{a}A{}^{a}\right)\Box \left(\nabla {}_{b}A{}^{b}\right)
%\end{align*}
%\section{Discussion}
%In \cite{} similar ideas and techniques will be applied to the graviton field in de Sitter.
%
%On the equivalence or not of the Landau and Lorentz gauge. We propose to use the Landau gauge even when for terms like
%$$
%\frac{z}{\lambda}
%$$ 
%Obviously the limits $\lambda \to \infty$ and $z to \infty$ do not commute in this situation. We propose however to apply the landau gauge in an ``agressive'' manner thus eliminating the term $\frac{z}{\lambda}$. When applied to the graviton field
\section{acknowledgments}
It is a pleasure to thank  J.P Gazeau, A. Higuchi, K. Meissner, U. Moschella, E. Mottola, A.M. Polyakov  and A. Starobinsky for stimulating discussions  about de Sitter quantum field theory. 
%\section{Appendix: Calculation of the flat space two-point function and propagator}
%Our conventions are 
%$$
%g_{aa'}=\text{diag}(-1,1,1,1), \quad f(x)=\int \frac{d^4k}{(2\pi)^4} e^{i k x} f(k)
%$$
%then
%the equation of motion reads
%$$
%\left[ g_{ab} \Box+(\lambda-1) \nabla_a \nabla_b\right] A^b(x)=0
%$$
%In Fourier space this is
%$$
%\left[ k^2 g_{ab} +(\lambda-1) k_a k_b\right] A^b(k)=0
%$$
%We define the Feynman propagator by
%\begin{align*}
%\langle \Omega\vert T \phi(x) \phi(x')\vert\Omega\rangle&=-i G^F(x,x')\\
%G^F(x,x')&=\int \frac{d^4k}{(2\pi)^4} \frac{e^{i k (x-x')}}{k^2+m^2-i \epsilon}\\
%(-\Box+m^2) G^F(x,x')&=\delta(x-x')
%\end{align*}
%We get 
%$$
%\left[ k^2 g_{a}^{b} +(\lambda-1) k_a k^b\right]  G^F_{b a'}=g_{aa'}
%$$
%with 
%$$
%G^F_{aa'}(x,x')=\frac{1}{k^2-i\epsilon} g_{aa'}+\frac{1-\lambda}{\lambda} \frac{1}{(k^2-i\epsilon)^2}k_a k_{a'}
%$$
%We have
%\begin{align*}
%\mathcal{G}_m(\mu)&=\int \frac{d^4k}{(2\pi)^4} \frac{e^{i k x}}{k^2+m^2-i \epsilon}=i\int_{\R^4} \frac{d^4k_E}{(2\pi)^4} \frac{e^{i k_E x}}{k_E^2+m^2}=\frac{i m}{(2\pi)^2} \frac{K_1(m \mu)}{\mu}\\
%\mathcal{G}_0(\mu)&=\frac{i}{4\pi^2 \mu^2}\\
%\mathcal{J}(\mu)&=\int \frac{d^4k}{(2\pi)^4} \frac{e^{i k x}}{(k^2-i \epsilon)^2}=i \int_{\R^4} \frac{d^4k_E}{(2\pi)^4} \frac{e^{i k_E x}}{k_E^4}=\frac{-i}{8 \pi^2} \ln \mu\\
%\end{align*}
%and
%\begin{align*}
%\omega_m(\mu)&=-i \mathcal{G}_m(\mu)=\frac{ m}{(2\pi)^2} \frac{K_1(m \mu)}{\mu}\\
%\omega_0(\mu)&=-i \mathcal{G}_0(\mu)=\frac{1}{4\pi^2 \mu^2}\\
%\varpi(\mu)&=-i \mathcal{J}(\mu)=\frac{-1}{8 \pi^2} \ln \mu
%\end{align*}
%
%Note that $\mathcal{J}$ is defined via it's derivative because the integral is not convergent.
%Using the formula
%$$
%\partial_a \partial_{a'} f(\mu)=\left( -\frac{f'}{\mu}\right) g_{aa'}+\left(-\frac{f'}{\mu}+f'' \right) n_{a} n_{a'}
%$$
%We find 
%\begin{align*}
%W_{aa'}&=\omega_0 g_{aa'}+\frac{1-\lambda}{\lambda} \varpi k_a k_{a'}\\
%&= \left[\frac{\lambda +1}{8 \pi ^2 \lambda  \mu ^2}\right] g_{aa'}-\left[\frac{\lambda -1}{4 \pi ^2 \lambda  \mu ^2}\right] n_{a} n_{a'}
%\end{align*}
%For the higher derivative theory we have
%$$
%\left[\Box g_{ab}+(\lambda-1) \nabla_a \nabla_b +\frac{\xi}{m^2} \Box  \nabla_a \nabla_b \right]A^b=0
%$$
%In Fourier space the equation verified by the Feynman propagator reads
%$$
%\left[k^2 g_{a}^{\:\:b}+(\lambda-1) k_a k^b-\frac{\xi}{m^2} k^2 k_a k^b \right] \tilde{G}^F_{ba'}=g_{a a'}
%$$  
%It is straightforward to invert this equation and one gets
%\begin{align*}
%\tilde{G}^F_{a a'}(k)=&\frac{ g_{aa'}}{k^2-i\epsilon}+\Bigg[\frac{1-\lambda }{\lambda\:  (k^2-i\epsilon)^2 } \\
%&+\frac{\xi}{m^2 \lambda ^2} \left( \frac{1}{k^2-i\epsilon}- \frac{1}{k^2-\frac{m^2 \lambda}{\xi}-i\epsilon}\right)\Bigg] k_a k_{a'}
%\end{align*}
%The reader will note the negative masse squared term appearing in the propagator which is typical of higher derivative theories. 
%We find 
%\begin{align*}
%W_{aa'}&= \left[ \frac{1}{2 \pi ^2 m^2 \lambda  \mu ^4 \nu ^2}+\frac{1}{2 \pi ^2 m^2 \lambda  \mu ^4 \nu ^2}+\frac{K_2(i m \mu  \nu )}{4 \pi ^2 \lambda  \mu ^2}\right] g_{aa'}+\\
%   &\left[\frac{2}{\pi ^2 m^2 \lambda  \mu ^4 \nu ^2}+\frac{1-\lambda }{4 \pi ^2 \lambda  \mu ^2}+\frac{K_0(i m \mu  \nu )}{\pi ^2 \lambda  \mu ^2}+\left(\frac{i m \nu }{4 \pi ^2 \lambda  \mu }-\frac{2 i}{\pi ^2 m
%   \lambda  \mu ^3 \nu }\right) K_1(i m \mu  \nu )\right] n_{a} n_{a'}
%\end{align*}
%The UV behavior is 
%$$
%W_{aa'}\approx \frac{1}{8\pi^2 \mu^2} \left[ g_{aa'}-2 n_a n_{a'}\right]
%$$
\begin{thebibliography}{breitestes Label}
\bibitem{Strocchi} F. Strocchi, \emph{Selected topics on the general properties of quantum field theory}, World Scientific Publishing, Singapore, $1993$.
\bibitem{MarolfHollands} A.M. Polyakov, \emph{De Sitter space and Eternity}, Nucl. Phys. B $797(2008)199$, A.M. Polyakov, \emph{Decay of vacuum energy}, Nucl. Phys. B $834(2010)316$, D. Marolf and I.A. Morrison, \emph{The IR stability of de Sitter: Loop corrections to scalar
propagators}, arXiv:1006.0035v3, D. Marolff and I.A. Morrison, \emph{The IR stability of de Sitter QFT: results at all orders} arXiv:1010.5327v1, S. Hollands, \emph{Correlators, Feynman diagrams, and quantum no-hair in
deSitter spacetime}, arXiv:1010.5367v1.
%\bibitem{Zuber} C. Itzykson and J-B. Zuber, \emph{Quantum field theory}, McGraw-Hill, Inc., New York, $1980$. 
\bibitem{Moschella} J. Bros, H. Epstein and U. Moschella, \emph{Scalar tachyons in the de Sitter universe}, Lett. math. phys. $93(2010)203$.
\bibitem{Staro79} A. A. Starobinsky, \emph{Spectrum of relict gravitational radiation and the early state of the universe}, JETP Lett. $30(1979)682$.
\bibitem{Allen86} B. Allen and T. Jacobson, \emph{Vector Two-Point Functions in Maximally Symmetric Spaces}, Commun. Math. Phys.  $103(1986)669$.
\bibitem{Antoniadis07} I. Antoniadis, P. Mazur and E. Mottola, \emph{Cosmological dark energy: prospects for a dynamical theory}, New J. Phys. $9(2007)11$.
\bibitem{Antoniadis91} I. Antoniadis and E. Mottola, \emph{Graviton fluctuations in de Sitter space}, J. Math. Phys. $32(1991)1037$
\bibitem{Higuchi} A. Higuchi and S.S. Kouris, \emph{The covariant graviton propagator in de Sitter spacetime}, 	Class. Quant. Grav. $18(2001)4317$.
%\bibitem{Woodard07}, S. P. Miao, N. C. Tsamis, and R. P. Woodard, \emph{de Sitter breaking through infrared divergences}, J. Math. Phys. $51(2010) 072503$.
\bibitem{ilio86} I Antoniadis, J. Iliopoulos and T.N. Tomaras, \emph{Gauge invariance in quantum gravity}, Nucl. Phys. B $267(1986)497$.
\bibitem{ilio87} E.G. Floratos, J. Iliopoulos and  T.N. Tomaras,\emph{Tree-level scatering amplitudes in de Sitter space diverge}, Phys. Lett. B $197(1987)373$
\bibitem{xAct} J.M. Martin-Garcia, \emph{xAct, Efficient Tensor Computer Algebra},  $2002$-$2008$
http://metric.iem.csic.es/Martin-Garcia/xAct/.
\bibitem{Ahmed} A. Youssef, \emph{Infrared behavior and gauge artifacts in the de Sitter universe II. The graviton field}, in preparation.
\end{thebibliography} 

\end{document}
