Table Of ContentNumerical approach to multi-loop integrals
2
1
0
2
n
a
J
K. Kato∗a†, E.deDonckerb,N.Hamaguchic, T.Ishikawac,T.Koiked, Y.Kuriharac,
0
3 Y. Shimizuc,andF.Yuasac
a KogakuinUniversity,1-24Nishi-Shinjuku,Shinjuku,Tokyo163-8677,Japan
]
h b WesternMichiganUniversityKalamazoo,MI49008-5371,USA
p cHighEnergyAcceleratorResearchOrganization(KEK),1-1OhoTsukuba,Ibaraki305-0801,
-
p Japan
e d SeikeiUniversity,Musashino,Tokyo180-8633,Japan
h
[
1
v Forthecalculationofmulti-loopFeynmanintegrals,anovelnumericalmethod,theDirectCom-
7
putationMethod(DCM)isdeveloped.Itisacombinationofanumericalintegrationandaseries
2
1 extrapolation. Inprinciple,DCM canhandlediagramsofarbitraryinternalmassesandexternal
6
momenta, and can calculate integrals with general numerator function. As an example of the
.
1
performanceofDCM,anumericalcomputationoftwo-loopboxdiagramsispresented. Further
0
2 discussion is given on the choice of controlparametersin DCM. Thismethodwill be an indis-
1
pensabletoolforthehigherorderradiativecorrectionwhenitistestedforawiderclassofphysical
:
v
parametersandtheseparationofdivergenceisdoneautomatically.
i
X
r
a
TheXXthInternationalWorkshopHighEnergyPhysicsandQuantumFieldTheory
September24-October1,2011
SochiRussia
∗Speaker.
†E-mail:[email protected]
PoS(QFTHEP2011)029
Numericalapproachtomulti-loopintegrals K.Kato
1. Introduction
Thehigh-statistics datainhigh-energy physicsrequiresthetheoreticalpredictionwithenough
accuracy. The prediction can be given by perturbative calculation in quantum field theory. Then
themulti-loop integralisanindispensable componentforthetheoretical study.
Wedefinethemulti-loopintegralsbythefollowingformulawherethespace-timedimensionis
denoted asn=4−2d 1. Hereweconfinethediscussion toscalarintegrals only. Sincethemethod
presented belowisbasicallynumerical, theinclusion ofanumeratorwillbestraight-forward.
I = (cid:213) L dnℓj (cid:213) N 1 (1.1)
Z (2p )ni D
j=1 r=1 r
wherethepropagator isD =q2−m2+ie ,N isthenumberofpropagators andListhenumberof
r r r
loops. Wecombinethepropagators bythestandardFeynmanparameterintegral
(cid:213) N 1 =(N−1)! (cid:213) dx d (1−(cid:229) xr) (1.2)
D Z r ((cid:229) x D )N
r=1 r r r
andperform integration withrespecttotheloopmomenta. Weobtain
I = G (N−nL/2)×I, I=(−1)N (cid:213) dx d (1−(cid:229) xr) , (1.3)
(4p )nL/2 Z rUn/2(V −ie )N−nL/2
V =M2−W, M2=(cid:229) x m2 (1.4)
U r r
r
whereU andW arepolynomials inthexparameters[1].
In Section 2, we propose a unique method to calculate the integral I in Eq.1.3. We call the
methodtheDirectComputationMethod(DCM).Intheprecedingworks[2],DCMhassuccessfully
calculated one-loop and two-loop diagrams. As an example we show the results for the two-loop
box diagrams in Section 3 and also report a study on the parameters in DCM. We discuss further
aspects ofDCMinSection4.
2. Method
DCM consists of a regularized numerical integration and an extrapolation of a numerical se-
quence.
In order to illustrate the idea of regularized numerical integration, let us consider a simple
integral:
1 dx
J= (2.1)
Z m2−sx(1−x)−ie
0
AsisshowninFig.1,theintegrand ofEq.2.1hassingular points whens>4m2. Analytically,
thiscanbehandledbytakinge asaninfinitesimalpositivequantity,or,inotherwords,bythehyper-
functionformula1/(z−ie )=P(1/z)+ipd (z). Numericalcomputationisunstableiftheintegrand
is divergent at some points. A way to solve the situation is to deform the path in the complex x
1Thesymbole isreservedforthe(infinitesimal)parameterinthepropagator.
2
Numericalapproachtomulti-loopintegrals K.Kato
x +ie
1 0
(a) (b) x +ie
1 1
x +ie
1 2
x x
0 1 0 1
Figure1: Integralpathsfor(a)analyticalcalculationandfor(b)DCM.Thecross(×)standsforthesingu-
larityofintegrand.
plane to avoid the singular points. However, the deformation would be very complicated for the
integrandofamulti-variablefunction. Anotherpossibilityistoassumethate isafinitequantityas
isshowninFig.1(b). Then,thenumericalintegration alongthex-axisisstable. Ifwecancalculate
thelimitingvaluefore →0,itisthevalueoftheintegralJ.
TheloopintegralI(e )isdefinedby
I(e )=(−1)N (cid:213) dx d (1−(cid:229) xr)(V +ie )N−nL/2, (2.2)
Z r Un/2(V2+e 2)N−nL/2
and if e and d are finite, the numerical integration can be performed using a suitable numerical
computation library.
Weusee determinedbya(scaled)geometric sequence
e =e =e /(A )l, (l=0,1,···) (2.3)
l 0 c
forconstants e ,A (A >1). Thenweexpect
0 c c
I= limI(e ). (2.4)
l
l→¥
Repeating the numerical integration, we obtain a sequence of numerical values of I(e ) for
l
l =0,1,...,l . Fromthese values andusing anextrapolation method, wecanestimate the value
max
ofI withenoughaccuracy.
Next we discuss the singularity originating from d →0. In Eq.2.2, V2+e 2 is positive and
U is positive semi-definite. Only at the boundary of the integration regionU becomes 0 2. When
1 1
d →+0,itiseitherintegrablelike dxdy ornon-integrablelike dx . Inthelatter
Z0 (x+y)1−d Z0 x1−d
case it develops a pole term ∼1/d as the ultraviolet singularity3. Depending on the masses and
external momenta, V can be 0 inside the integration region to develop the imaginary part of I 4,
and also can be 0 at the boundary of integral region (as in the case ofU) to develop an infrared
2U isasumofmonomialsofx.
3ThesingularpolealsoappearsinthefirstfactorofEq.1.3ifN−nL/2≤0ford =0.
4Forillustration,oneassumesthatN−nL/2=1andthevariablesaretransformedintoV andx′ variables. Then,
r
omittingtheJacobianandotherdetails,theimaginarypartbecomes (cid:213) dx′ V2dV e . Inthelimitase →+0,
Z rZV1 V2+e 2
theinnerintegralisfiniteifV2>0>V1and0otherwise.
3
Numericalapproachtomulti-loopintegrals K.Kato
singularity pole ∼1/d . If Eq.2.2 is free from these singularities, we just put d =0. If not, the
integral in Eq.2.2 is denoted as I(d ,e ) and we first calculate I(d )=liml→¥ I(d ,el) as in Eq.2.4
fixingd . Then,weassumethefollowingform:
C
I(d )=···+ −1 +C +C d +··· (2.5)
d 0 1
Wecalculate I(d )forasetofvalues ofd andestimate thecoefficientsC . Forinstance, incaseof
j
asinglepole,d I(d )=C−1+C0d +O(d 2)andweextractC−1 andC0 5.
So much for the description of DCM and one can understand the necessity of an efficient
libraryforthenumericalintegrationandthatfortheextrapolation. FortheformerweuseDQAGE[3]
whichisavariant ofGaussian quadrature. Sinceitworksadaptively, onecanspecify theaccuracy
of the numerical results, although the high accuracy costs in computation time. For the latter we
use Wynn’s e algorithm[4] 6 which predicts the limiting value by the following iteration. We set
theresultsofthenumericalintegration asinitialvaluesoftheseriesa(l,k):
a(l,−1)=0, a(l,0)=I(e ), l=0,1,···. (2.6)
l
Theelementa(l,k+1)isobtainedbythefollowingrecurrence relation:
1
a(l,k+1)=a(l+1,k−1)+ , l=0,1,···. (2.7)
a(l+1,k)−a(l,k)
Whilst the a(l,k) values with odd k are meant to store temporary numbers, the a(l,k) with even k
giveextrapolated estimates.
3. Numerical results
Wecalculate the integrals for the two-loop box diagrams shown in Fig.2. Theparameters are
m =m =m =m =m=50GeV, m =m =m =M=90GeV, p2 = p2 = p3 = p2 =m2 and
1 2 5 6 3 4 7 1 2 3 4
t =(p +p )2 =−(100)2GeV2. We take s=(p +p )2 variable and introduce a dimensionless
1 3 1 2
variable f =s/m2.
s
(a) (b)
p1 1 5 p3 p1 1 5 p3
3 4 7 3 4 7
p2 2 6 p4 p2 2 6 p4
Figure2: (a)Two-loopplanarboxdiagramand(b)Two-loopnon-planarboxdiagram.Themassofinternal
linekism . Eachexternalmomentumflowsinward.
k
TheexplicitformoftheU andW functionsisfoundin[5]. Sincethereisnoultraviolet/infrared
divergence,weputd =0inEq.1.3. Theintegralis6-dimensionalandweperformatransformation
5AndalsoC isnecessaryifthefirstfactorofEq.1.3issingular.
1
6This’e ’hasnothingtodowiththeparameterinpropagators.
4
Numericalapproachtomulti-loopintegrals K.Kato
oftheintegration variablesontoa6-dimensional hypercube [0,1]6. Bythistransformation onecan
cancelcommonvariablesbetweenthenumeratorandthedenominator. Theresultsarepresentedin
Fig.3andinFig.4. In[5],weverifiedtheresultsbycomparisonwithanothercomputationwhichis
acombination ofalgebraic transformations andnumerical integration, andalsobytheconsistency
checkbetweentherealandimaginarypartsthrough thedispersion relation.
3
Re
Im
2.5
2
1.5
) 1
s,t
(
I 0.5
0
-0.5
-1
-1.5
0 5 10 15 20 25
f
s
Figure 3: Numericalresultsof I forthe planardiagramin unitsof10−12 GeV−6 for0.0≤ f ≤25.0and
s
t=−10000.0GeV2. Plottedpointsaretherealpart(bullets)andtheimaginarypart(squares).Forthelatter,
theregion f >10isnotyetcomputed.
s
Re
0.8 Im
0.6
0.4
)
s,t
(
I
0.2
0
-0.2
0 5 10 15 20
f
s
Figure4: Numericalresultsof I forthenon-planardiagramin unitsof 10−12 GeV−6 for0.0≤ f ≤20.0
s
andt=−10000.0GeV2. Plottedpointsaretherealpart(bullets)andtheimaginarypart(squares).
In DCM, the validity of the extrapolation depends on the choice of the values of e . We keep
finite value of e in the numerical integration, so that its physical dimension is the same as the
squared mass. Wehavetwoparameters inEq.2.3. A isnormally A =2andwecan useA =1.2
c c c
or1.3toobtainalesscomputational intensivesequenceofintegralsase decreasesmoreslowly. In
order to study the choice of e , the initial value of the iteration, we have calculated the following
0
5
Numericalapproachtomulti-loopintegrals K.Kato
5
4
] 3
2
−
V
e 2
G
5
−
0 1
1
)[ M=90GeV
eI( 0 Realpart
Imaginarypart
-1
-2
0.01 1 100 10000 1e+06
e
Figure 5: The real part(lower points) and imaginary part(upperpoints) of one-loop scalar vertex integral
I(e )areshown.Thehorizontalredlineistheexact(analytical)value.
on-mass-shell one-loop vertexintegralasanexample:
1
I= dxdy (3.1)
Z0≤x+y≤1 M2(1−x−y)+m2e(x+y)2−sxy−ie
Here,s=5002GeV2andm =0.5×10−3GeV. Thevalueofthisintegraliscomputedwithagiven
e
valueofM ande =1.270−m form=0,1,...,120. Thenweusethevaluesform,m+1,...,m+14
as the input of extrapolation. This means that we take l =14 and e =1.270−mGeV2 for m=
max 0
0,1,...,100.
Table.1 Extrapolated values of I. 15 values, I(e )= e /(1.2)l, (l = 0,...,14), are used for the
l 0
extrapolation. Errorisnotthedifference fromtheanalytical valuebutestimated fromtheextrapo-
lation.
e RealpartofI error ImaginarypartofI error
0
3.49E+05 -1.75104242540072E-05 1.65E-09 2.25556360020931E-09 1.80E-09
5.63E+04 -1.75105247961553E-05 3.02E-13 -2.96789310051170E-05 8.78E-08
9.10E+03 -1.75105248407057E-05 1.55E-14 4.35402513918612E-05 1.82E-09
1.47E+03 -1.75105250993412E-05 6.78E-21 4.34982757936633E-05 4.04E-13
2.37E+02 -1.75105250996915E-05 2.76E-16 4.34982757936633E-05 4.04E-13
3.83E+01 -1.75105250989527E-05 1.40E-17 4.34982788205288E-05 1.92E-17
6.19E+00 -1.75105250987944E-05 3.58E-18 4.34982788206820E-05 5.55E-16
1.00E+00 -1.75105251104117E-05 1.22E-15 4.34982788654878E-05 9.93E-14
1.62E-01 -1.75105250921691E-05 1.98E-16 4.34982793676347E-05 2.11E-14
2.61E-02 -1.75105251352866E-05 1.71E-15 4.34982790401757E-05 4.97E-14
4.21E-03 -1.75105251110498E-05 4.07E-15 4.34982788582295E-05 2.21E-17
(analytical) -1.75105250974494E-05 4.34982788194091E-05
Thecalculated results forM=90GeVareshowninFig.5andTable.1. Itistobenoted thate
inDCMisobviouslyfinite. OnecanseethatevenwhenthevalueofI(e )differsfromtheanalytical
value, the extrapolation gives agood estimation. Thereisafinite region thatshowsthe agreement
between the analytical value and the extrapolated one. This behavior demonstrates that DCM is
stable up to some extent for the choice of e parameters. We have tested the similar analysis for
0
6
Numericalapproachtomulti-loopintegrals K.Kato
several values ofM (M =1,101,102,103GeV)and found similarbehavior. Though weneedmore
testsforthispoint,wecantemporallyconcludethate isnotneedtobeverysmallbutitcanbeset
0
tothetypicalsquaredmassinthedenominator oftheintegrand.
Inthecalculationofthetwo-loopboxdiagramsdescribedabove,wehavecheckedtheconver-
gence ofthe extrapolation step-by-step. The value used for e is1.240 ∼1.245GeV2 which would
0
beconsistent withtheaboveconjecture.
4. Summary
In this paper, we have outlined DCM and calculated the two-loop scalar box integral as an
exampletoshowitsapplicability. Sincetheradiativecorrectionintheelectroweaktheory(orinthe
SUSYmodel)involves various combinations ofmassparameters intheintegrand, DCMisagood
candidate tohandlegeneralloopintegrals.
Inorder touse DCMfor thecalculation ofhigher-order radiative corrections, weplan topro-
ceedtothefollowingresearch.
1. The method should be tested for a wider class of diagrams with various combinations of
massesandexternalmomentaandwiththenumerator structure.
2. Furtherstudyonthechoiceofparameters e ,A isrequiredforastableapplication.
0 c
3. The variable transformation is sometimes important for good convergence. This is to be
processed inanautomaticmanner.
4. Indimensional regularization, theultraviolet/infrared divergence appearsasapole1/d . The
separation of the infrared pole is already done successfully in [6]. A similar treatment of
ultraviolet polesisexpected.
5. Sometimes DCM needs long computational time. It will be important to perform the com-
putations inaparallelcomputing environment.
Acknowledgements
WewishtothankProf.T.Kanekoforhisdiscussionsandvaluablecomments.Thisworkwassupported
inpartbytheGrant-in-Aid(No.20340063andNo.23540328)ofJSPSandbytheCPISprogramofSokendai.
References
[1] R.J.Eden,P.V.Landshoff,D.I.Olive,J.C.Polkinghorne,TheAnalyticS-Matrix,CambridgeUniversity
Press,Cambridge,1966;
G.Tiktopoulos,Phys.RevVol131(1963)480.
[2] See[5]andreferencestherein.
[3] R.Piessens,E.deDoncker,C.W.UbelhuberandD.K.Kahaner,QUADPACK,ASubroutinePackage
forAutomaticIntegration,SpringerSeriesinComputationalMathematics.Springer-Verlag,(1983).
[4] D.Shanks,J.Math.Phys.34(1955)1;
P.Wynn,MathematicalTablesAidstoComputing10(1956)91;SIAMJ.Numer.Anal.3(1966)91.
[5] F.Yuasa,E.deDoncker,N.Hamaguchi,T.Ishikawa,K.Kato,Y.Kurihara,J.Fujimoto,Y.Shimizu,
[arXiv:1112.0637/hep-ph],submittedtoComput.Phys.Commun.(2011).
[6] E.deDoncker,J.Fujimoto,N.Hamaguchi,T.Ishikawa,Y.Kurihara,M.Ljucovic,Y.Shimizu,F.Yuasa,
PoS(CPP2010)011[arXiv:1110.3587/hep-ph],(2011).
7
