restart;with(plots):## simplest form
beta:=unapply(-beta0*as^2,as);## initial condition
inc:=as(mz2)=alpha0;
## the (ordinary) differential equation for alpha_s
ODE:=Q2*diff(as(Q2),Q2)=beta(as(Q2));
## solving the ODE
as1:=dsolve({ODE,inc},as(Q2)) assuming Q2>0,mz2>0;
## doing some algebraic simplifications...
as2:=rhs(subs(Q2=Q^2,mz2=mz^2,combine(as1,ln))) assuming Q2>0,mz2>0;## define the Lambda parameter...
Lambda:=unapply(Q^2/exp(1/(beta0*as(Q^2))),Q);
## ... express alpha0 (the initial condition from above) in terms of Lambda
t1:=solve(subs(as(mz^2)=a0,Lambda(mz))=Lambda^2,a0) assuming Q>0,mz>0;
## and express alphas in terms of Lambda
as013:=1/combine(simplify(subs(alpha0=t1,1/as2)),ln) assuming Q>0,mz>0,Lambda>0;## Plot the solution for different values for Lambda
## Note: There is a factor of 4 Pi to convert as -> alphas
beta0:=11-2/3*nf;
nf:=5;
p1:=semilogplot(subs(Lambda=0.350,4*Pi*as013),Q=10..2000,axes=boxed,color=red,labels=["Q [GeV]",alphas],labeldirections=[horizontal, vertical],thickness=1,view=[10..2000,0.05..0.25],legend="Lambda = 350 MeV",tickmarks=[21,10],numpoints=1000):
p2:=semilogplot(subs(Lambda=0.200,4*Pi*as013),Q=10..2000,axes=boxed,color=blue,labels=["Q [GeV]",alphas],labeldirections=[horizontal, vertical],thickness=1,legend="Lambda = 200 MeV",numpoints=1000):
p3:=semilogplot(subs(Lambda=0.100,4*Pi*as013),Q=10..2000,axes=boxed,color=green,labels=["Q [GeV]",alphas],labeldirections=[horizontal, vertical],thickness=1,legend="Lambda = 100 MeV",numpoints=1000):
p4:=semilogplot(subs(Lambda=0.050,4*Pi*as013),Q=10..2000,axes=boxed,color=magenta,labels=["Q [GeV]",alphas],labeldirections=[horizontal, vertical],thickness=1,legend="Lambda = 50 MeV",numpoints=1000):
display({p1,p2,p3,p4});
## now higher orders....
unassign('beta0','as');
beta:=-(beta0*as^2+beta1*as^3+beta2*as^4+beta3*as^5);
##2-loop:
## expand the 1/beta into a Taylor series...
as10:=convert(series(1/beta,as,4),polynom) assuming as>0;
## integrate ...
as11:=(int(as10,as)+C)=L;
## substitute C = beta1/beta0^2*log(beta0) and do some simplifiactions...
as111:=expand(beta0*combine(subs(C=beta1/beta0^2*log(beta0),as11),ln)) assuming L>0,beta0>0;
## solve for the 1/as - term and substitute the 1-loop solution into the other as-dependent terms
## and expand into a 1/L - series
as112:=series(subs(L=1/l,subs(as = 1/beta0/L,solve(subs(1/as=1/a,as111),a))),l,3);
as113:=combine(subs(l=1/L,convert(as112,polynom)),ln) assuming L>0,beta0>0;##3-loop:
as20:=convert(series(1/beta,as,5),polynom) assuming as>0:
as21:=(int(as20,as)+C)=L:
as211:=expand(beta0*combine(subs(C=beta1/beta0^2*log(beta0),as21),ln)) assuming L>0,beta0>0:
as212:=series(subs(L=1/l,subs(as = as113,solve(subs(1/as=1/a,as211),a))),l,4):
as213:=combine(subs(l=1/L,convert(as212,polynom)),ln) assuming L>0,beta0>0:
as214:=collect(as213,L);##4-loop:
as30:=convert(series(1/beta,as,6),polynom) assuming as>0:
as31:=(int(as30,as)+C)=L:
as311:=expand(beta0*combine(subs(C=beta1/beta0^2*log(beta0),as31),ln)) assuming L>0,beta0>0:
as312:=expand(series(subs(L=1/l,subs(as = as214,(solve(subs(1/as=1/a,as311),a)))),l,5)):
as313:=simplify(subs(l=1/L,convert(as312,polynom)),ln) assuming L>0,beta0>0:
as314:=collect(as313,L);## some 4-loop plots...
beta0:=11-2/3*nf;
beta1:=102-38*nf/3;
beta2:=(77139-15099*nf+325*nf^2)/54;
beta3:=29243-6946.3*nf+405.089*nf^2+1.49931*nf^3;L:=log(Q^2/Lambda^2):
## range for the plots...
xmin:=10:
xmax:=2000:
ymin:=0.06:
ymax:=0.22:
p1:=semilogplot(subs(Lambda=0.350,nf=5,4*Pi*as314),Q=10..2000,axes=boxed,color=red,labels=["Q [GeV]",alphas],labeldirections=[horizontal, vertical],thickness=1,view=[xmin..xmax,ymin..ymax],legend="Lambda = 350 MeV",tickmarks=[21,10],numpoints=1000):
p2:=semilogplot(subs(Lambda=0.200,nf=5,4*Pi*as314),Q=10..2000,axes=boxed,color=blue,labels=["Q [GeV]",alphas],labeldirections=[horizontal, vertical],thickness=1,legend="Lambda = 200 MeV",numpoints=1000):
p3:=semilogplot(subs(Lambda=0.100,nf=5,4*Pi*as314),Q=10..2000,axes=boxed,color=green,labels=["Q [GeV]",alphas],labeldirections=[horizontal, vertical],thickness=1,legend="Lambda = 100 MeV",numpoints=1000):
p4:=semilogplot(subs(Lambda=0.050,nf=5,4*Pi*as314),Q=10..2000,axes=boxed,color=magenta,labels=["Q [GeV]",alphas],labeldirections=[horizontal, vertical],thickness=1,legend="Lambda = 50 MeV",numpoints=1000):
display({p1,p2,p3,p4});## Choose Lambda such that - order by order - as(mz) = 0.1134 (0.1184)
unassign('Lambda');
## the Z^0 mass
mz:=91.18;
## initial value alpha_s(Q=mz), needed for the following computations,
## change it and reexecute the worksheet
a0:=0.1134;
Lambda_0:=abs(fsolve(subs(Q=mz,nf=5,4*Pi*as013) = a0,Lambda));
Lambda_1:=abs(fsolve(subs(Q=mz,nf=5,4*Pi*as113) = a0,Lambda));
Lambda_2:=abs(fsolve(subs(Q=mz,nf=5,4*Pi*as214) = a0,Lambda));
Lambda_3:=abs(fsolve(subs(Q=mz,nf=5,4*Pi*as314) = a0,Lambda));
## Fractional differences... 1 - as(n-loop)/as(4-loop)
nf:=5;
## range for the plots...
xmin:=10:
xmax:=2000:
ymin:=-0.002:
ymax:=0.002:
f1:=semilogplot(1 - subs(Lambda=Lambda_0,as013)/subs(Lambda = Lambda_3,as314),Q=10..2000,axes=boxed,color=red,labels=["Q [GeV]",""],labeldirections=[horizontal, vertical],thickness=1,view=[xmin..xmax,ymin..ymax],legend="1-loop",tickmarks=[21,10],numpoints=1000):
f2:=semilogplot(1 - subs(Lambda=Lambda_1,as113)/subs(Lambda = Lambda_3,as314),Q=10..2000,axes=boxed,color=blue,labels=["Q [GeV]",""],labeldirections=[horizontal, vertical],thickness=1,view=[xmin..xmax,ymin..ymax],legend="2-loop",tickmarks=[21,10],numpoints=1000):
f3:=semilogplot(1 - subs(Lambda=Lambda_2,as214)/subs(Lambda = Lambda_3,as314),Q=10..2000,axes=boxed,color=green,labels=["Q [GeV]",""],labeldirections=[horizontal, vertical],thickness=1,view=[xmin..xmax,ymin..ymax],legend="3-loop",tickmarks=[21,10],numpoints=1000):
f4:=semilogplot(0,Q=10..2000,axes=boxed,color=black,labels=["Q [GeV]",""],labeldirections=[horizontal, vertical],thickness=0,view=[xmin..xmax,ymin..ymax],legend="4-loop",tickmarks=[21,10],numpoints=1000):
display({f1,f2,f3,f4});## now with quark threshold matching...
## matching condition at three loop: as(nf-1)/as(nf) = 1 + C2 as(nf)^2 + C3 as(nf) at the pole mass
## first step...
unassign('nf');
C2:=-0.291667;
##C3:=-5.32389+(nf-1)*0.26247;
C3:=0;
a1:= (a/Pi)*(1 + C2*(a/Pi)^2 +C3*(a/Pi)^3);
beta0:=11-2/3*nf;
beta1:=102-38*nf/3;
beta2:=(77139-15099*nf+325*nf^2)/54;
beta3:=29243-6946.3*nf+405.089*nf^2+1.49931*nf^3;
## the heavy quark masses
Mbottom:=4.53;
Mcharm:=1.5;
Lambda5:=Lambda_3;## calculate as at the bottom threshold and match it to nf = 4
a5:=evalf(subs(Q=Mbottom,nf=5,Lambda=Lambda5,4*Pi*as214));
a15:=evalf(Pi*subs(a=a5,nf=5,a1));
a5/a15;## calculate Lambda such that as(nf=4, Q = Mbottom) = a15
Lambda4:=abs(fsolve(subs(Q=Mbottom,nf=4,4*Pi*as214) = a15,Lambda));
#evalf(subs(Q=Mbottom,nf=4,Lambda=0.234,4*Pi*as214));
#Lambda4:=0.2955;## calculate as at the charm threshold and match it to nf = 3
a4:=evalf(subs(Q=Mcharm,nf=4,Lambda=Lambda4,4*Pi*as214));
a14:=evalf(Pi*subs(a=a4,nf=4,a1));
a4/a14;#evalf(subs(Q=Mcharm,nf=3,Lambda=0.3385,4*Pi*as214));
#Lambda3:=0.3385;
Lambda3:=abs(fsolve(subs(Q=Mcharm,nf=3,4*Pi*as214) = a14,Lambda));## some plots...
h1:=plot(subs(Lambda=Lambda3,nf=3,4*Pi*as214),Q=1.2..Mcharm,color=red,axes=boxed,labels=["Q [GeV]",alphas],labeldirections=[horizontal, vertical],thickness=1,view=[1.2..Mcharm,0.3..0.35],numpoints=1000):
h2:=plot(subs(Lambda=Lambda4,nf=4,4*Pi*as214),Q=Mcharm..Mbottom,color=blue,axes=boxed,labels=["Q [GeV]",alphas],labeldirections=[horizontal, vertical],thickness=1,view=[Mcharm..Mbottom,0.2..0.35],numpoints=1000):
h3:=plot(subs(Lambda=Lambda5,nf=5,4*Pi*as214),Q=Mbottom..10,color=green,axes=boxed,labels=["Q [GeV]",alphas],labeldirections=[horizontal, vertical],thickness=1,view=[Mbottom..7,0.15..0.2],numpoints=1000):
h:=plot(subs(Lambda=Lambda5,nf=5,4*Pi*as314),Q=1.2..10,color=black,axes=boxed,labels=["Q [GeV]",alphas],labeldirections=[horizontal, vertical],thickness=0,view=[Mbottom..7,0.15..0.2],numpoints=1000):
display({h1,h2,h3,h});j1:=plot(subs(Lambda=Lambda3,nf=3,4*Pi*as214)/subs(Lambda=Lambda5,nf=5,4*Pi*as314)-1,Q=1.1..Mcharm,color=red,axes=boxed,labels=["Q [GeV]",""],labeldirections=[horizontal, vertical],thickness=1,view=[1.1..Mcharm,-0.01..0.12],numpoints=1000):
j2:=plot(subs(Lambda=Lambda4,nf=4,4*Pi*as214)/subs(Lambda=Lambda5,nf=5,4*Pi*as314)-1,Q=Mcharm..Mbottom,color=blue,axes=boxed,labels=["Q [GeV]",""],labeldirections=[horizontal, vertical],thickness=1,view=[Mcharm..Mbottom,-0.01..0.12],numpoints=1000):
j3:=plot(subs(Lambda=Lambda5,nf=5,4*Pi*as214)/subs(Lambda=Lambda5,nf=5,4*Pi*as314)-1,Q=Mbottom..10,color=green,axes=boxed,labels=["Q [GeV]",""],labeldirections=[horizontal, vertical],thickness=1,view=[Mbottom..7,-0.01..0.12],numpoints=1000):
display({j1,j2,j3});