We start by loading FeynCalc and the necessary extensions

$LoadAddOns={"FeynArts", "FeynOnium","FeynHelpers"};
<<FeynCalc`
$FAVerbose = 0;

Configuration

We will be using the t'Hooft-Veltman-Breitenlohner-Maison scheme to deal with (γ^5) in D  dimensions. In this case no finite counter-term is needed.

FCSetDiracGammaScheme["BMHV"];

For pedagogical reasons will explicitly distinguish between UV and IR divergences at 1-loop

$KeepLogDivergentScalelessIntegrals=True;

This just makes the typesetting of the momenta look better

MapThread[FCAttachTypesettingRule[#1,{SubscriptBox,"p",#2}]&,{{p1,p2},Range[2]}];
MapThread[FCAttachTypesettingRule[#1,{SubscriptBox,"k",#2}]&,{{k1,k2},Range[2]}];

Generate Feynman diagrams

We use FeynArts to generate the tree-level and 1-loop diagrams for the process c (Overscript[c, -]) -> γ γ

diagsTree = InsertFields[
CreateTopologies[0,2->2],(*2->2 process*)
{F[3,{2,SUNFIndex[cq]}],-F[3,{2,SUNFIndex[cqbar]}]}->{V[1],V[1]}, (*ccbar to gamma gamma*)
InsertionLevel->{Particles},
Model->"SMQCD", (*QCD extension of the EW Standard Model*)
ExcludeParticles->{S[_],V[1|2|3|4],U[_],F[1|2|4],F[3,{1|3}]} (*Remove irrelevant particles (Higgs, W,Z, other quark flavors etc.) in the internal lines *)
];

Visualize the generated diagrams

Paint[diagsTree, ColumnsXRows -> {2, 1}, Numbering -> None,SheetHeader->None,ImageSize->{512,256}];

Same story for the 1-loop diagrams

diagsLoop = InsertFields[
CreateTopologies[1,2->2,ExcludeTopologies->{Tadpoles,WFCorrections}],{F[3,{2,cq}],-F[3,{2,cqbar}]}->{V[1],V[1]},
InsertionLevel->{Particles},
Model->"SMQCD",
ExcludeParticles->{S[_],V[1|2|3|4],U[_],F[1|2|4],F[3,{1|3}]}
];
Paint[DiagramExtract[diagsLoop,3..10],ColumnsXRows->{4,2}, Numbering->Simple,
SheetHeader->False,ImageSize->{1024,512}];

Get the amplitudes

CreateFeynAmp is a FeynArts function that generates the amplitude
FCFAConvert is a FeynCalc function that converts this amplitude to the FeynCalc notation

ampTree[0]=FCFAConvert[
CreateFeynAmp[diagsTree,Truncated->False,PreFactor->1],
IncomingMomenta->p,
OutgoingMomenta->k,
UndoChiralSplittings->True,(*no need to split each vertex into left and right chiral parts*)
TransversePolarizationVectors->{k1,k2},
ChangeDimension->D,(*generate D-dimensional amplitude*)
SMP->True,
FinalSubstitutions->{SMP["m_c"]->SMP["m_Q"]},
Contract->True];

Replace the electric charge of the quark by (Subscript[e, Q]), introduce the fine-structure constant α

ampTree[1]=Factor[(3/2 SMP["e_Q"])^2 ampTree[0]/.SMP["e"]^2->SMP["alpha_fs"]4Pi]

ampLoop[0]=FCFAConvert[CreateFeynAmp[DiagramExtract[diagsLoop,3..10],
Truncated->False,PreFactor->1],
IncomingMomenta->p,OutgoingMomenta->k,LoopMomenta->{l},
UndoChiralSplittings->True,
TransversePolarizationVectors->{k1,k2},
ChangeDimension->D,
SMP->True,
DropSumOver->True,(*remove summations over external colors as introduced by FeynArts and used in FormCalc*)
FinalSubstitutions->{SMP["m_c"]->SMP["m_Q"]},Contract->True];

ampLoop[1]=Factor[(3/2 SMP["e_Q"])^2 ampLoop[0]/.{SMP["e"]^2->SMP["alpha_fs"]4Pi,
SMP["g_s"]^2->SMP["alpha_s"]4Pi}]

Set up the kinematics

Reset the existing kinematic  constraints

FCClearScalarProducts[]

Just a convenient shortcut

QMass=SMP["m_Q"];

Putting all external particles on the mass shell

ScalarProduct[p1,p1]=QMass^2;
ScalarProduct[p2,p2]=QMass^2;
ScalarProduct[k1,k1]=0;
ScalarProduct[k2,k2]=0;

Additional rest frame constraints. We also introduce "hat"-vectors here

ScalarProduct[P,q]=0;
ScalarProduct[P,P]=4 EN^2;
ScalarProduct[k1,P]=2EN^2;
ScalarProduct[k2,P]=2EN^2;
ScalarProduct[k2,q]=-FCI[SP[k1,q]];
ScalarProduct[k1,k2]=2 EN^2;

ScalarProduct[k1hat,k1hat]=0;
ScalarProduct[qhat,qhat]=-1;
ScalarProduct[Phat,q]=0;
ScalarProduct[Phat,Phat]=1;
ScalarProduct[k1hat,P]=2EN;
ScalarProduct[k1hat,Phat]=1;
ScalarProduct[k2hat,Phat]=1;
ScalarProduct[Phat,qhat]=0;

ScalarProduct[P,Polarization[k2,-I,Transversality->True]]=0;
ScalarProduct[P,Polarization[k1,-I,Transversality->True]]=0;

Finally, we can also specify temporal and spatial components of the relevant vectors

TC[q]=0;
TC[P]=2EN;
TemporalMomentum[Polarization[k1|k2,___]]=0;
CartesianMomentum[P|Phat,___]=0;
CartesianMomentum[k2hat,dim___]=-CartesianMomentum[k1hat,dim];
CSP[qhat,qhat]=1;
CSP[k1hat,k1hat]=1;
TC[k1]=EN;
TC[k1hat]=1;
TC[k2hat]=1;
TC[Phat]=1;
TC[qhat]=0;

CartesianScalarProduct[k1hat,Polarization[k1,-I,Transversality->True]]=0;
CartesianScalarProduct[k1hat,Polarization[k2,-I,Transversality->True]]=0;

This is a small hack to avoid simplifying expressions like (Sqrt[(m^2)]) by hand

Unprotect[Power];
Power[Power[QMass,2],Rational[1,2]]:=QMass;
Protect[Power];

Tree-level amplitudes

First of all, we flag the charm mass in the propagator to distinguish it from (Subscript[m, c]) that comes from putting external momenta on-shell.
This is irrelevant for a pure tree-level calculation but quite important for a loop calculation, where we need to renormalize (Subscript[m, c]).
Only the charm mass from the propagator must be renormalized!

The second command inserts the covariant-projector that we need to extract the contribution with the proper quantum numbers. We can immediately
take the Dirac trace thereafter.

ampTree[2]=ampTree[1]/.QMass->mark[QMass]/.Spinor[x__]:>(Spinor[x]/.mark->Identity)//
FMInsertCovariantProjector[#,{p1,QMass},{p2,QMass},{cq,cqbar},Dimension->D]&//SUNSimplify//DiracSimplify

Here we switch to the rest-frame 2-body kinematics with P and q.

ampTree[3]=ampTree[2]//FCReplaceMomenta[#,{p1->1/2P+q,p2->1/2P-q,k2->P-k1}]&//
	ReplaceAll[#,x_FeynAmpDenominator/;FreeQ[{x},l]:>FeynAmpDenominatorExplicit[x]]&//
	ExpandScalarProduct[#,EpsEvaluate->True]&;

Replace P by 2*E*(OverHat[P])  and (Subscript[k, 1]) by E*(OverHat[(Subscript[k, 1])]) . Also rewrite  E as  and q as |(Overscript[q, -> ])|*(OverHat[q])

ampTree[4]=ampTree[3]//ReplaceRepeated[#,{Momentum[P,dim___]->2EN Momentum[Phat,dim],
Momentum[k1,dim___]:>EN Momentum[k1hat,dim],EN->Sqrt[QMass^2-FCI@SP[q,q]],
	Momentum[q,dim___]:>qabs Momentum[qhat]}]&;

Expand in |(Overscript[q, -> ])| up to second order and do the replacement (q^μ)(q^ν)-> (((Overscript[q, -> ])^2)/3)(-(g^(μ ν))+ (((P^μ)(P^ν))/(P^2)))

ampTree[5]=Series[ampTree[4],{qabs,0,2}]//Normal//Uncontract[#,qhat,Pair->All]&//
ReplaceAll[#,FCI[FV[qhat,i_]FV[qhat,j_]]:>FCI[CSP[qhat,qhat]/3*(-MT[i,j]+FV[P,i]FV[P,j]/SP[P,P])]]&//
Contract//ReplaceAll[#,Momentum[P]->2EN Momentum[Phat]]&;

Add up both amplitudes and switch to the Cartesian form

(Total[ampTree[5]]/.mark->Identity)//LorentzToCartesian

One-loop amplitudes

Tensor decomposition to the basis of Passarino-Veltman coefficient functions

AbsoluteTiming[
counter=0;
ampLoop[2]=(counter++; WriteString["stdout","Amp ",counter, " out of ", Length[ampLoop[1]]," ...\n"];TID[#,l,ToPaVe->True,UsePaVeBasis->True])&/@ampLoop[1];
]

Insert the covariant projector, calculate the traces

AbsoluteTiming[
counter=0;
ampLoop[3]=(counter++; WriteString["stdout","Amp ",counter, " out of ", Length[ampLoop[1]]," ...\n"];SUNSimplify[DiracSimplify[FMInsertCovariantProjector[#,{p1,QMass},{p2,QMass},{cq,cqbar},Dimension->D]]])&/@ampLoop[2];
]

Switch to the rest-frame 2-body kinematics with P and q.

AbsoluteTiming[ampLoop[4]=ampLoop[3]//FCReplaceMomenta[#,{p1->1/2P+q,p2->1/2P-q,k2->P-k1}]&//
	ReplaceAll[#,x_FeynAmpDenominator/;FreeQ[{x},l]:>FeynAmpDenominatorExplicit[x]]&//
	ExpandScalarProduct[#,EpsEvaluate->True]&;]

Replace P by 2*E*(OverHat[P])  and (Subscript[k, 1]) by E*(OverHat[(Subscript[k, 1])]) . Also rewrite  E as  and q as |(Overscript[q, -> ])|*(OverHat[q])

AbsoluteTiming[ampLoop[5]=ampLoop[4]//ReplaceRepeated[#,{Momentum[P,dim___]->2EN Momentum[Phat,dim],
Momentum[k1,dim___]:>EN Momentum[k1hat,dim],EN->Sqrt[QMass^2-FCI@SP[q,q]],
	Momentum[q,dim___]:>qabs Momentum[qhat]}]&;]

Using Package-X we can automatically expand the Passarino-Veltman functions in |(Overscript[q, -> ])| up to the desired order and evaluate them
analytically

AbsoluteTiming[
counter=0;
ampLoop[6]=(counter++; WriteString["stdout","Amp ",counter, " out of ", Length[ampLoop[1]]," ...\n"];PaXEvaluateUVIRSplit[#,l,
PaXImplicitPrefactor->1/(2Pi)^D,PaXSeries->{{qabs,0,2}},PaXAnalytic->True])&/@ampLoop[5];
]

Do the replacement (q^μ)(q^ν)-> (((Overscript[q, -> ])^2)/3)(-(g^(μ ν))+ (((P^μ)(P^ν))/(P^2)))

AbsoluteTiming[ampLoop[7]=Uncontract[ampLoop[6],qhat,Pair->All]//Collect2[#,qhat]&//ReplaceAll[#,FCI[FV[qhat,i_]FV[qhat,j_]]:>FCI[CSP[qhat,qhat]/3](-MT[i,j]+FV[P,i]FV[P,j]/SP[P,P])]&//Contract//
ReplaceAll[#,Momentum[P]->2EN Momentum[Phat]]&//FCHideEpsilon//SimplifyPolyLog//ReplaceAll[#,Zeta2->Zeta[2]]&;]

One last missing ingredient is the UV-renormalization. To that aim we multiply the tree-level amplitude by the wave-function renormalization
constant in the on-shell scheme and renormalize the charm mass in the same scheme

ruleRenormalization={
Zpsi->1-CF SMP["alpha_s"]/(4Pi)(SMP["Delta_UV"]+2SMP["Delta_IR"]+3Log[ScaleMu^2/QMass^2]+4),
	mark[QMass]->(1-3*CF*SMP["alpha_s"]*(4/3 + SMP["Delta_UV"]+Log[ScaleMu^2/QMass^2])/(4*Pi))QMass
}

ampLoopWFCorrection=(Total[ampTree[5]](Zpsi)/.ruleRenormalization)//FCShowEpsilon//Series[#,{SMP["alpha_s"],0,1}]&//Normal//FCHideEpsilon

Now we are essentially done. Dropping the IR-divergences that will cancel against those in the NRQCD amplitude we have

ampFull=((Total[ampLoop[7]]+ ampLoopWFCorrection)/.SMP["Delta_IR"]->0)//ChangeDimension[#,4]&

It might seem like the amplitude has a contribution proportional to a single power of v.

Coefficient[ampFull,qabs,1]

However, upon rewriting the 4D-LeviCivita tensor in terms of Cartesian 3D epsilon tensors, we see that this contribution vanishes

Coefficient[ampFull,qabs,1]//LorentzToCartesian//ExpandAll

So let us extract the values of c0 and c2

overallPrefactor=(Sqrt[2]Sqrt[CA] FCI@LC[][k1hat,Phat,Polarization[k1,-I,Transversality->True],Polarization[k2,-I,Transversality->True]])

c0Coeff=Collect2[Cancel[Coefficient[ampFull,qabs,0]/overallPrefactor],SMP["alpha_s"]]//Cancel

c2Coeff=Collect2[Cancel[QMass^2 Coefficient[ampFull,qabs,2]/(overallPrefactor)],SMP["alpha_s"]]

And so we arrive to the final matching coefficients

matchingCoeffF=Collect2[QMass^2/(8 Pi)c0Coeff^2//Series[#,{SMP["alpha_s"],0,1}]&//Normal,SMP["alpha_s"],FCFactorOut->2Pi SMP["e_Q"]^4 SMP["alpha_fs"]^2]

matchingCoeffG=Collect2[QMass^2/(4 Pi)(c0Coeff c2Coeff)//Series[#,{SMP["alpha_s"],0,1}]&//Normal,SMP["alpha_s"],FCFactorOut->-8Pi/3 SMP["e_Q"]^4 SMP["alpha_fs"]^2,Factoring->ExpandAll]

Compare this to the results obtained in arXiv:1104.1418 Eq. 38

matchingCoeffF$Literature=2 Pi SMP["e_Q"]^4 SMP["alpha_fs"]^2(1+ CF SMP["alpha_s"]/Pi(Pi^2/4-5));
matchingCoeffG$Literature=-8/3 Pi SMP["e_Q"]^4 SMP["alpha_fs"]^2(1+ CF SMP["alpha_s"]/Pi(5Pi^2/16-49/12-Log[ScaleMu^2/(4 QMass^2)]));

FCCompareResults[matchingCoeffF,matchingCoeffF$Literature]
FCCompareResults[matchingCoeffG//PowerExpand,matchingCoeffG$Literature//PowerExpand]