One-loop gravity with xAct

In this notebook, we will do the explicit computation of the one-loop diagrams contributing to the two-point function in gravity. In the following, we will (mostly) ignore factors of i. We will not compute the contribution from the Faddeev-Popov ghosts, which are needed to get a diffeomorphism-invariant result.

Reading the package

Block[{Print},
<<xAct`xTras`;
];

Defining a manifold and metric perturbations

Most of this should look familiar. New elements are:
1) the definition of momenta and their squares, plus automatic replacement rules (the notation is [thing you want to assign the rule to]/:[rule])
2) some “nice” constant symbols to compute the propagator (you can check what a command does by evaluating the command with a question mark in front: ?DefNiceConstantSymbol)

Block[{Print},
DefManifold[spacetime,4,{μ,ν,ρ,σ,α,β,γ,δ,κ,λ,τ,ω}];
DefMetric[-1,metric[-μ,-ν],CD,{";","𝒟"}];
DefMetricPerturbation[metric,pert,ϵ,PrintAs->"h"];
DefTensor[hpert[-μ,-ν],spacetime,Symmetric[{-μ,-ν}],PrintAs->ColorString["𝒽",Lighter@Blue]];
metricpertrule={pert[LI[1],X__]:>hpert[X],pert[LI[n_/;n>1],__]:>0};

DefConstantSymbol[GNewton,PrintAs->"(Subscript[G, N])"];

DefConstantSymbol[GFα,PrintAs->"(Subscript[α, gf])"];
DefConstantSymbol[GFβ,PrintAs->"(Subscript[β, gf])"];

DefTensor[p1[-μ],spacetime,PrintAs->"(Subscript[𝓅, 1])"];
p1/:CD[μ_]@p1[ν_]:=0;
DefTensor[loopmom[-μ],spacetime,PrintAs->"ℓ"];
loopmom/:CD[μ_]@loopmom[ν_]:=0;

DefConstantSymbol[p1sq,PrintAs->"[(Subscript[𝓅, 1]^2)]"];
DefConstantSymbol[loopmomsq,PrintAs->"[(ℓ^2)]"];
DefConstantSymbol[p1plusloopmomsq,PrintAs->"[((Subscript[𝓅, 1])+ℓ)^2]"]; 

p1/:p1[-μ_]p1[μ_]:=p1sq;
p1/:p1[μ_]loopmom[ν_]/;(μ===-ν):=((p1plusloopmomsq-p1sq-loopmomsq)/2);
loopmom/:loopmom[-μ_]loopmom[μ_]:=loopmomsq;

Table[DefNiceConstantSymbol[G,i],{i,1,5}];
];

?DefNiceConstantSymbol

You might notice that we defined the scalar product of p1 and loopmom in a somewhat involved way. The reason is that denominators are tricky to simplify, and using (((p+q))^2) as a scalar is much better than expanding the expression, since this is what appears in denominators.

Computing the propagator

Let us define a general gauge fixing with an arbitrary parameter β:

ℱgf[α_,μ_,ν_]:=((1/2)(metric[α,μ]CD[ν]@#+metric[α,ν]CD[μ]@#)-((1+GFβ)/4)metric[μ,ν]CD[α]@#)&;

In a general gauge, the graviton propagator can have five different terms. This is due to symmetry: it must be symmetric in each of the pairs, and also symmetric under exchange (convince yourself of that!). The overall behaviour of the propagator should be 1/(p^2), which we take into account in our ansatz.

proph[μ_,ν_,ρ_,σ_,pmom_]:=Module[{momsq=CollectTensors[pmom[α]pmom[-α]]},(1/momsq)((1/2)G1(metric[μ,ρ]metric[ν,σ]+metric[μ,σ]metric[ν,ρ])+G2 metric[μ,ν]metric[ρ,σ]+(G3/4)(1/momsq)(metric[μ,ρ]pmom[ν]pmom[σ]+metric[μ,σ]pmom[ν]pmom[ρ]+metric[ν,ρ]pmom[μ]pmom[σ]+metric[ν,σ]pmom[μ]pmom[ρ])+(G4/2)(1/momsq)(metric[μ,ν]pmom[ρ]pmom[σ]+metric[ρ,σ]pmom[μ]pmom[ν])+(G5/(momsq^2)) pmom[μ]pmom[ν]pmom[ρ]pmom[σ])];

The first step to compute the propagator is to compute the quadratic part of the action, and to add the gauge fixing. We will introduce a second gauge parameter α, indicating how strongly we impose the gauge. The limit α->0 is called the Landau limit, and corresponds to a strict implementation of the gauge condition.

S2=(1/(16π GNewton))CollectTensors[ToFlat[ExpandPerturbation[((Perturbation[(Sqrt[-Detmetric[]])RicciScalarCD[],2])/(Sqrt[-Detmetric[]]))]/.metricpertrule]-(1/(GFα))ℱgf[α,μ,ν]@hpert[-μ,-ν]ℱgf[-α,ρ,σ]@hpert[-ρ,-σ]];

We next take one functional derivative so that all covariant derivatives act on one of the fields:

dS2dh=VarD[hpert[κ,λ],CD][S2]//CollectTensors;

We will compute diagrams in a flat background, so we can go to momentum space. This entails replacing covariant derivatives by the momentum:

dS2dhmom=dS2dh//.CD[μ_]@hpert[X__]:>I p1[μ]hpert[X]//CollectTensors;

What is this expression? It is the two-point function acting on the field h. If we now multiply the propagator from the left, we should get back the field. This is the defining property of the propagator - it is the inverse of the two-point function.

propdS2dhmom=proph[α,β,κ,λ,p1]dS2dhmom//CollectTensors;

We can use that this expression should simply be the field to compute the coefficients of the propagator:

prophcoeffs=SolveConstants[propdS2dhmom==hpert[α,β],{G1,G2,G3,G4,G5}][[1]]//Simplify

Now we have the propagator. To compute loop diagrams, we still need vertices. We will compute both the self-energy diagram and the tadpole, even though the latter vanishes in dimreg.
The two-point diagrams will have two pairs of external indices, call them (μν) and (ρσ). We also have a single external momentum, call it p1.
We will also have internal indices, which we label (αβ), (γδ), (κλ), (τω). The loop momentum will be called loopmom.

Computing the 3h vertices

For the self-energy diagram, we need two 3h-vertices. They come with exactly the opposite signs of momenta.
How do we compute these vertices? In particular, how can we transform covariant derivatives to momenta in this case (since we have three different momenta)? For this, we will code our own version of VarD, called VD, that works directly in Fourier space. Let’s have a bunch of definitions first:

VD[tenlist_?ListQ][momlist_?ListQ,alist_?ListQ][exprlist_?ListQ]:=Module[{counting},Table[VD[tenlist][momlist,alist][exprlist[[counting]]],{counting,1,Length@exprlist}]];
VD[tenlist_?ListQ][momlist_?ListQ,alist_?ListQ][expr_]:=If[Length@momlist==1,VD[tenlist[[1]]][momlist[[1]],alist[[1]]][expr],VD[Drop[tenlist,1]][Drop[momlist,1],Drop[alist,1]][CollectTensors@VD[tenlist[[1]]][momlist[[1]],alist[[1]]][expr]]];
VD[ten_?xTensorQ][mom_,a_][CD[μ_]@X__]:=I mom[μ]VD[ten][mom,a][X];
VD[ten_?xTensorQ][mom_,a_][expr_/;Head[expr]===Plus]:=Plus@@Table[VD[ten][mom,a][expr[[Q]]],{Q,1,Length@expr}];
VD[ten_?xTensorQ][mom_,a_][expr_/;Head[expr]===Times]:=Plus@@Table[VD[ten][mom,a]@expr[[Q]] Drop[expr,{Q}],{Q,1,Length@expr}];
VD[ten_?xTensorQ][mom_,a_][expr_?ConstantQ]=0;
VD[ten_?xTensorQ][mom_,a_][delta[_,_]]=0;
VD[ten_?xTensorQ][mom_,a_][metric[__]]=0;
VD[ten_?xTensorQ][mom_,a_][p1[__]]=0;
VD[ten_?xTensorQ][mom_,a_][loopmom[__]]=0;
VD[hpert][mom_,{μ_,ν_}][hpert[ρ_,σ_]]:=(1/2)(metric[-μ,ρ]metric[-ν,σ]+metric[-μ,σ]metric[-ν,ρ]);

Let’s explain this mess. First of all, what is the syntax of VD? VD needs three angular brackets: the first is either a single tensor, or a list of tensors, with respect to which we take the functional derivative(s). VD is written in such a way that it can take multiple derivatives in one go. The second set of arguments is a momentum and an index (or a list of both, if we take multiple derivatives; arguments have to be in the same order as the list of tensors). Lastly, the third angular brackets enclose the thing that we take the derivative(s) of.

So far, so good. Now for the actual implementation:
- The first line implements what is called “Listability”: if the final argument is a list, VD simply applies to each list element, and we get back a list.
- The second line implements the multiple-derivatives-at-once-part: it is a recursive definition that repeats VD until we take a functional derivative w.r.t. a single field in each application.
- The third line implements the Fourier transform: covariant derivatives get transformed into momenta.
- The fourth line implements that acting with VD on a sum gives the sum of VD acting on all the summands.
- The fifth line implements the product rule.
- The sixth line implements that the functional derivative of a constant expression is zero.
- The lines 7 to 10 implement that functional derivatives of the Dirac delta, the metric, and the momenta p1 and loopmom vanish.
- Finally, the last line implements that taking the functional derivative of h w.r.t. h is the identity (on the space of symmetric rank two tensors).

We can now compute the actual vertices from the expansion of the action to cubic order in h:

CollectTensors[(1/(16π GNewton))ToFlat[ExpandPerturbation[((Perturbation[(Sqrt[-Detmetric[]])RicciScalarCD[],3])/(Sqrt[-Detmetric[]]))]/.metricpertrule]];
vert1=VD[{hpert,hpert,hpert}][{p1,loopmom,(-p1[#]-loopmom[#])&},{{μ,ν},{α,β},{κ,λ}}][%]//CollectTensors;
vert2=VD[{hpert,hpert,hpert}][{(-p1[#])&,(-loopmom[#])&,(p1[#]+loopmom[#])&},{{ρ,σ},{γ,δ},{τ,ω}}][%%]//CollectTensors;

Now, what the he...misphere are the momentum arguments doing? The key issue is that in our definition of VD, we assume that the arguments are pure tensors. This means that if we take this argument and throw a bracket with an index on it, like p1[μ], this has to be a valid expression. Now, for some of the momentum arguments, we either need a minus sign for momentum conservation, or even the sum of two momenta. The problem is now that (p1+loopmom)[μ] is not a valid expression (it is simply not implemented in this way).
This is where the “#&” construction comes into play - this is called a “pure function”. Basically, the point is that a pure function can act on any argument (that is, putting angular brackets with arguments in it after the pure function). The result is that the argument is substituted where the hash is in the pure function. This means e.g. that the pure function (#+1)& (which you could call “add one”), acting on e.g. 5, gives 6: (#+1)&[5] -> 6.
This construction works with our definition of VD: the argument we supply is an index, and it simply gets substituted as the index of all the momenta. Problem solved!

Computing the self-energy diagram

We now have all ingredients to compute the self-energy diagram. The only thing that is still left is to multiply everything and to contract indices. It is helpful to go step by step (this speeds up the computation).
We first contract one vertex with one propagator:

V1G=vert1 proph[α,β,γ,δ,loopmom]/.prophcoeffs//CollectTensors;

Next, we multiply the second vertex with the other propagator:

V2G=vert2 proph[κ,λ,τ,ω,(loopmom[#]+p1[#])&]/.prophcoeffs//CollectTensors;

Let us now compute the product of these two expressions. It is a bit faster to split this product into a list of products, and then combine terms again (it still takes a while!).

Monitor[selfenergytab=Table[V1G[[i]]V2G//CollectTensors,{i,1,Length@V1G}];,Row[{i,"/",Length@V1G}]]

selfenergy=CollectTensors[Total[selfenergytab]];

As said before, the final expression will carry four indices. We however know that this expression has the same symmetries as the propagator, so we can span it in terms of the same basis. The formulas to do so are derived in the solutions to sheet 3.
For the implementation of this rule, it is important to go from many to fewer loop momenta, since the rule is applied in order. This means we first replace the terms with four loop momenta, then three, and so on. Think about why this is the right way!
We can use MakeRule to simplify our lives a bit. This... makes a replacement rule. Amazing, right? Note that on the RHS of the rule, we need to use Evaluate, because otherwise Mathematica “holds” the expression (in brief: it doesn’t work otherwise).

MakeRule[{loopmom[μ]loopmom[ν]loopmom[ρ]loopmom[σ],(1/5)(loopmomsq^2)(((1-(((([(ℓ^2)]-[((((Subscript[𝓅, 1])+ℓ))^2)]+[(Subscript[𝓅, 1]^2)]))^2)/(4 [(ℓ^2)] [(Subscript[𝓅, 1]^2)]))))^2)Evaluate[CollectTensors[Symmetrize[metric[μ,ν]metric[ρ,σ]]]]-(2/5)((loopmomsq^2)/p1sq)(1-7(((([(ℓ^2)]-[((((Subscript[𝓅, 1])+ℓ))^2)]+[(Subscript[𝓅, 1]^2)]))^2)/(4 [(ℓ^2)] [(Subscript[𝓅, 1]^2)]))+6(((([(ℓ^2)]-[((((Subscript[𝓅, 1])+ℓ))^2)]+[(Subscript[𝓅, 1]^2)]))^4)/(16 (([(ℓ^2)])^2) (([(Subscript[𝓅, 1]^2)])^2))))Evaluate[CollectTensors[Symmetrize[metric[μ,ν]p1[ρ]p1[σ]]]]+(1/5)((loopmomsq^2)/(p1sq^2))(1-12(((([(ℓ^2)]-[((((Subscript[𝓅, 1])+ℓ))^2)]+[(Subscript[𝓅, 1]^2)]))^2)/(4 [(ℓ^2)] [(Subscript[𝓅, 1]^2)]))+16(((([(ℓ^2)]-[((((Subscript[𝓅, 1])+ℓ))^2)]+[(Subscript[𝓅, 1]^2)]))^4)/(16 (([(ℓ^2)])^2) (([(Subscript[𝓅, 1]^2)])^2))))p1[μ]p1[ν]p1[ρ]p1[σ]}];
MakeRule[{loopmom[μ]loopmom[ν]loopmom[ρ],(loopmomsq/p1sq)((p1plusloopmomsq-p1sq-loopmomsq)/2)(1-(((([(ℓ^2)]-[((((Subscript[𝓅, 1])+ℓ))^2)]+[(Subscript[𝓅, 1]^2)]))^2)/(4 [(ℓ^2)] [(Subscript[𝓅, 1]^2)])))Evaluate[CollectTensors[Symmetrize[metric[μ,ν]p1[ρ]]]]-(loopmomsq/(p1sq^2))((p1plusloopmomsq-p1sq-loopmomsq)/2)(1-2(((([(ℓ^2)]-[((((Subscript[𝓅, 1])+ℓ))^2)]+[(Subscript[𝓅, 1]^2)]))^2)/(4 [(ℓ^2)] [(Subscript[𝓅, 1]^2)])))p1[μ]p1[ν]p1[ρ]}];
MakeRule[{loopmom[μ]loopmom[ν],(1/3)loopmomsq(1-(((([(ℓ^2)]-[((((Subscript[𝓅, 1])+ℓ))^2)]+[(Subscript[𝓅, 1]^2)]))^2)/(4 [(ℓ^2)] [(Subscript[𝓅, 1]^2)])))metric[μ,ν]-(1/3)(loopmomsq/p1sq)(1-4(((([(ℓ^2)]-[((((Subscript[𝓅, 1])+ℓ))^2)]+[(Subscript[𝓅, 1]^2)]))^2)/(4 [(ℓ^2)] [(Subscript[𝓅, 1]^2)])))p1[μ]p1[ν]}];
MakeRule[{loopmom[μ],(1/p1sq)((p1plusloopmomsq-p1sq-loopmomsq)/2) p1[μ]}];
sloppyloopmomrules=Join[%%%%,%%%,%%,%];

Let us apply this rule:

selfenergyfinal=CollectTensors[selfenergy/.sloppyloopmomrules];

The result is still a complete mess of course, even in a specific gauge:

selfenergyfinal/.GFα->1/.GFβ->1

Computing the 4h vertices

Let us also compute the tadpole diagram. For this, we need the four-graviton vertex. In complete analogy:

CollectTensors[(1/(16π GNewton))ToFlat[ExpandPerturbation[((Perturbation[(Sqrt[-Detmetric[]])RicciScalarCD[],4])/(Sqrt[-Detmetric[]]))]/.metricpertrule]];
vert4h=VD[{hpert,hpert,hpert,hpert}][{p1,(-p1[#])&,loopmom,(-loopmom[#])&},{{μ,ν},{ρ,σ},{α,β},{γ,δ}}][%]//CollectTensors;

The full diagram contracts the vertex with a single propagator (this takes some time again):

tadpole=CollectTensors[vert4h proph[α,β,γ,δ,loopmom]/.prophcoeffs];

We once again apply our sloppy rule:

tadpolefinal=CollectTensors[tadpole/.sloppyloopmomrules];

What is missing

To get to the final result, there are a few steps left. First, we have to compute the contribution of the Faddeev-Popov ghosts. We then add up all the diagrams with their corresponding symmetry factors, and perform the loop integrals (remember that we didn’t take care of factors of i, and you would also have to do a Wick rotation). The resulting logarithmic divergence should be proportional to (p^4), and the tensor structure should be the same as the one of the two-point function of quadratic curvature terms.