Please note that before starting you need to fix the Working precision (wpc). Normaly 30 is enough;

Same for AccuracyGoal (acc = 5 is enough).

wpc=40;
$MinPrecision=wpc;
acc=7;

Entanglement entropies of disks in confining backgrounds

General equations

Equation that we want to solve

The equation that we want to solve is:

equation=ρ''[R]+((((ρ'[R])^3)(R^3))/2)((2Q-4)f[R] + R f'[R])+(4-Q)(((ρ'[R])^2)/(ρ[R])) + ((Q/R)+((f'[R])/(f[R])))ρ'[R]+((4-Q)/((R^4)f[R]ρ[R]))

Let us express q in terms of the confining scale (by solving f((Subscript[R, conf])) = 0.

Solve[1-(μ/(R^(Q-1)))-((q^2)/(R^(2Q-4)))==0,q]/.R->Rconf//Simplify

Using ζ = (Subscript[R, conf])/R, now f becomes

rulef={f->Function[ζ,Assuming[Q>0&&Rt>0&&ζ>0,1-(μ/(R^(Q-1)))-((q^2)/(R^(2Q-4)))/.q->(Rconf^(-2+(Q/2))) (Sqrt[(Rconf^Q)-Rconf μ])/.μ->μt (Rconf^(Q-1))/.R->Rconf/ζ//Simplify]//Evaluate]}

Let us now work out the change of variables between R and ζ:

dζdR = (D[Rconf/R,R]//Simplify)/.R->Rconf/ζ//Simplify;
changeRtoζ={func_''[R]-> dζdR D[dζdR func'[ζ],ζ],func_'[R]-> dζdR func'[ζ],func_[R]->func[ζ]}//Simplify

Then the dimensionless equation is:

equation/.changeRtoζ/.ρ-> ((Rconf^(-1)) ρt[#]&)/.R->Rconf/ζ//Simplify
Collect[(Rconf^3)(ζ^(-4))%,{ρt[ζ],ρt''[ζ],ρt'[ζ]},FullSimplify]

We factor out the dimensionful quantities:

equation/.changeRtoζ/.ρ-> ((Rconf^(-1)) ρt[#]&)/.R->Rconf/ζ/.rulef//Simplify;
Solve[%==0,ρt''[ζ]]//Flatten//Simplify;
equationdimless=%[[1,1]]-%[[1,2]]

Some parameters for the numerical routine

For numerics, it is convenient to define the equation in terms of its parameters:

Clear[eqNum]
eqNum[Q_,μt_]=equationdimless;

ϵIRnow=(10^(-5));ϵUVnow=(10^(-5));
plotrange={{-1,1},{0,-1}};
aspectRatio=1;
imageSize=300;
labelStyle=Directive[Black,FontFamily-> "Arial",FontSize->14];

We define the following function to plot the solutions:

plotSol[sol_]:=
(If[sol[[0]]=!= List,Print["solution not found"];Return[];];
ParametricPlot[{{-ρt[ζ]/.sol,-ζ}
,{ρt[ζ]/.sol,-ζ}},{ζ,sol[[1,2,1,1,1]],sol[[1,2,1,1,2]]},Ticks->{{{-1,-1},{1,1}},{{0,0},{-1,1}}},PlotRange->plotrange,AspectRatio->aspectRatio,ImageSize->imageSize,PlotStyle->Thick,LabelStyle->labelStyle]
)

Boundary conditions: Embedding for small spheres

orcusp1=1;

solcusp1={cR[0]->((σ1 2  (Sqrt[-3+Q]) ζs)/(Sqrt[2 (-2+Q) ζs-(-3+Q) ((ζs)^Q) μt])),cR[1]->-(((-2+Q) (2 (5+(-4+Q) Q) ((ζs)^5)+(((-3+Q))^3) ((ζs)^(4+Q)) μt-(((-3+Q))^2) (-2+Q) ((ζs)^(3+2 Q)) ((μt)^2)-2 (-3+Q) (-1+Q) ((ζs)^(1+2 Q)) (-1+μt) ((σ1)^2)+(((-3+Q))^2) (-1+Q) ((ζs)^(3 Q)) (-1+μt) μt ((σ1)^2)))/(2 (Sqrt[-3+Q]) (-1+Q) (((2 (-2+Q) ζs-(-3+Q) ((ζs)^Q) μt))^(3/2)) (((ζs)^4)+((ζs)^(2 Q)) (-1+μt)-((ζs)^(3+Q)) μt) σ1))};

Here σ1 is either + or minus; it orresponds to just picking one of the branches. It is enough to pick one of them, due to the symmetry with respect to ρ = 0. The plus sign picks the ρt positive branch

orcusp1=1;

Assuming[0<ζ<ζs<1&&Q>2&&μt>0,(equationdimless/.{ρt->Function[ζ,(∑_(kk=0)^orcusp1)cR[kk](((ζs-ζ))^(kk+(1/2)))+((O[ζ,ζs])^(orcusp1+(3/2)))]}/.solcusp1//Simplify)/.((σ1)^(n_)):>If[EvenQ[n],1,σ1]//Simplify]
Solve[%==0,cR[orcusp1]];
Assuming[0<ζ<ζs<1&&Q>2&&-2 (-2+Q) ζs+(-3+Q) ((ζs)^Q) μt<0,%//Simplify]

One could in principle find more terms but it is actually not necessary. It just takes some time (due to not having specified Q I think.

Clear[ρtcust1,dρtcust1];

ρtcust1[Q_,μt_,ζs_,ζ_]=(∑_(kk=0)^orcusp1)cR[kk](((ζs-ζ))^(kk+(1/2)))/.solcusp1/.σ1->1;
dρtcust1[Q_,μt_,ζs_,ζ_]=D[%,ζ];

The following function will compute the solution and store it in solnumerical, for each choice of the parametes

computeSol[Q_,μt_,ζs_]:=(solnumerical[Q,μt,ζs]=NDSolve[{eqNum[Q,μt]==0,ρt[ζs-ϵIRnow]==ρtcust1[Q,μt,ζs,ζs-ϵIRnow],ρt'[ζs-ϵIRnow]==dρtcust1[Q,μt,ζs,ζs-ϵIRnow]},ρt,{ζ,ζs-ϵIRnow,ϵUVnow},WorkingPrecision->wpc,Method->"StiffnessSwitching"]//Flatten);

Here we are actually going way fancier than we need. We could just use big numbers in the boundary conditions, see for example

computeSol[5,0,1/2];

Block[{Q=5,μt=0,inf=(10^8),ζs=1/2,aux1},
aux1=NDSolve[{eqNum[Q,μt]==0,ρt[ζs-ϵIRnow]==(inf^(-1)),ρt'[ζs-ϵIRnow]==inf},ρt,{ζ,ζs-ϵIRnow,ϵUVnow},WorkingPrecision->wpc,Method->"StiffnessSwitching"]//Flatten;

Plot[(ρt[ζ]/.aux1)-ρt[ζ]/.solnumerical[5,0,1/2],{ζ,ζs-ϵIRnow,ϵUVnow},PlotRange->All]

]

So the error is order ~(10^(-3)). Note however that it gets bigger near the cusp, where ρt is actually small, meaning the absolute error is not so small.

The plot looks like:

computeSol[5,0,1/2]

plotSol[solnumerical[5,0,1/2]]

Boundary conditions: Embedding for two spheres (not considered in the paper)

solcusp2={};

orcusp2=0;

solcusp2={cR[1]->((σ2 2 ζs)/(Sqrt[2 (-2+Q) ζs-(-3+Q) ((ζs)^Q) μt])),cR[2]->((2 (-4+Q) ((ζs)^2))/(6 (-2+Q) ζs ρs-3 (-3+Q) ((ζs)^Q) μt ρs)),cR[3]->((-6 (-2+Q) ζs+6 (-3+Q) ((ζs)^Q) μt-3 (-3+Q) Q ((ζs)^Q) μt-(6 (2 (-2+Q) ζs-(-3+Q) ((ζs)^Q) μt) ((-2+Q) ((ζs)^4)-(-2+Q) ((ζs)^(2 Q)) (-1+μt)+((ζs)^(3+Q)) μt))/(((ζs)^4)+((ζs)^(2 Q)) (-1+μt)-((ζs)^(3+Q)) μt)-((4 (((-4+Q))^2) ((ζs)^3))/((ρs)^2))+((4 (-4+Q) (-11+2 Q) ((ζs)^3))/((ρs)^2))) σ2)/(18 (((2 (-2+Q) ζs-(-3+Q) ((ζs)^Q) μt))^(3/2))),cR[4]->((-4+Q) ζs (2 (298-77 Q+4 (Q^2)) ((ζs)^7)-2 (298-77 Q+4 (Q^2)) ((ζs)^(6+Q)) μt+18 (22-19 Q+4 (Q^2)) ((ζs)^5) ((ρs)^2)+18 (14-11 Q+2 (Q^2)) ((ζs)^(1+2 Q)) (-1+μt) ((ρs)^2)-9 (82-67 Q+13 (Q^2)) ((ζs)^(4+Q)) μt ((ρs)^2)-45 (6-5 Q+(Q^2)) ((ζs)^(3 Q)) (-1+μt) μt ((ρs)^2)+((ζs)^(3+2 Q)) (-596+596 μt+351 ((μt)^2) ((ρs)^2)+Q (154-154 μt-279 ((μt)^2) ((ρs)^2))+(Q^2) (-8+8 μt+54 ((μt)^2) ((ρs)^2)))))/(135 (((-2 (-2+Q) ζs+(-3+Q) ((ζs)^Q) μt))^2) (((ζs)^4)+((ζs)^(2 Q)) (-1+μt)-((ζs)^(3+Q)) μt) ((ρs)^3)),cR[5]->((16 (7504-4600 Q+945 (Q^2)-70 (Q^3)+(Q^4)) ((ζs)^14)-32 (7504-4600 Q+945 (Q^2)-70 (Q^3)+(Q^4)) ((ζs)^(13+Q)) μt+16 (7504-4600 Q+945 (Q^2)-70 (Q^3)+(Q^4)) ((ζs)^(2 (6+Q))) ((μt)^2)-48 (-888+882 Q-257 (Q^2)+15 (Q^3)+2 (Q^4)) ((ζs)^12) ((ρs)^2)-24 (5880-6410 Q+2291 (Q^2)-292 (Q^3)+7 (Q^4)) ((ζs)^(11+Q)) μt ((ρs)^2)+36 (((-2+Q))^2) (33-28 Q+4 (Q^2)) ((ζs)^10) ((ρs)^4)-72 (((-2+Q))^2) (279-272 Q+68 (Q^2)) ((ζs)^(6+2 Q)) (-1+μt) ((ρs)^4)+36 (((-2+Q))^2) (-15-4 Q+4 (Q^2)) ((ζs)^(2+4 Q)) (((-1+μt))^2) ((ρs)^4)+36 (-60-200 Q+417 (Q^2)-235 (Q^3)+42 (Q^4)) ((ζs)^(9+Q)) μt ((ρs)^4)+72 (2376-4212 Q+2786 (Q^2)-813 (Q^3)+88 (Q^4)) ((ζs)^(5+3 Q)) (-1+μt) μt ((ρs)^4)-36 (-12-128 Q+173 (Q^2)-73 (Q^3)+10 (Q^4)) ((ζs)^(1+5 Q)) (((-1+μt))^2) μt ((ρs)^4)+45 (((-3+Q))^2) (4-12 Q+5 (Q^2)) ((ζs)^(6 Q)) (((-1+μt))^2) ((μt)^2) ((ρs)^4)+90 (6-5 Q+(Q^2)) ((ζs)^(3+5 Q)) (-1+μt) μt ((ρs)^2) (Q (-44+44 μt-19 ((μt)^2) ((ρs)^2))+(Q^2) (4-4 μt+((μt)^2) ((ρs)^2))+16 (7-7 μt+3 ((μt)^2) ((ρs)^2)))-8 (-4+Q) ((ζs)^(9+3 Q)) μt ((Q^3) (-4+4 μt+39 ((μt)^2) ((ρs)^2))-6 (Q^2) (-44+44 μt+81 ((μt)^2) ((ρs)^2))-16 (-469+469 μt+108 ((μt)^2) ((ρs)^2))+3 Q (-908+908 μt+561 ((μt)^2) ((ρs)^2)))+16 (-4+Q) ((ζs)^(2 (5+Q))) (3752-3752 μt-2403 ((μt)^2) ((ρs)^2)+2 (Q^3) (-1+μt+18 ((μt)^2) ((ρs)^2))-6 (Q^2) (-22+22 μt+95 ((μt)^2) ((ρs)^2))+3 Q (-454+454 μt+733 ((μt)^2) ((ρs)^2)))+6 (-2+Q) ((ζs)^(7+3 Q)) μt ((ρs)^2) (-288 (61-61 μt+2 ((μt)^2) ((ρs)^2))+3 (Q^3) (40-40 μt+33 ((μt)^2) ((ρs)^2))-8 (Q^2) (280-280 μt+87 ((μt)^2) ((ρs)^2))+Q (11432-11432 μt+1389 ((μt)^2) ((ρs)^2)))-3 (-2+Q) ((ζs)^(8+2 Q)) ((ρs)^2) (645 (Q^3) ((μt)^2) ((ρs)^2)-32 (Q^2) (37-37 μt+105 ((μt)^2) ((ρs)^2))-6 (2624-2624 μt+219 ((μt)^2) ((ρs)^2))+Q (8672-8672 μt+4677 ((μt)^2) ((ρs)^2)))-6 (-2+Q) ((ζs)^(4+4 Q)) (-1+μt) ((ρs)^2) ((Q^3) (16-16 μt+381 ((μt)^2) ((ρs)^2))-8 (Q^2) (55-55 μt+432 ((μt)^2) ((ρs)^2))-18 (240-240 μt+547 ((μt)^2) ((ρs)^2))+Q (2584-2584 μt+10221 ((μt)^2) ((ρs)^2)))+((ζs)^(6+4 Q)) ((Q^4) (16-32 μt+672 ((μt)^3) ((ρs)^2)+9 ((μt)^4) ((ρs)^4)+((μt)^2) (16-672 ((ρs)^2)))-16 (-7504+15008 μt-7236 ((μt)^3) ((ρs)^2)+81 ((μt)^4) ((ρs)^4)+268 ((μt)^2) (-28+27 ((ρs)^2)))-9 (Q^2) (-1680+3360 μt-6784 ((μt)^3) ((ρs)^2)+151 ((μt)^4) ((ρs)^4)+16 ((μt)^2) (-105+424 ((ρs)^2)))+2 (Q^3) (-560+1120 μt-5448 ((μt)^3) ((ρs)^2)+81 ((μt)^4) ((ρs)^4)+8 ((μt)^2) (-70+681 ((ρs)^2)))+8 Q (-9200+18400 μt-17730 ((μt)^3) ((ρs)^2)+351 ((μt)^4) ((ρs)^4)+10 ((μt)^2) (-920+1773 ((ρs)^2))))) σ2)/(2160 ζs (((2 (-2+Q) ζs-(-3+Q) ((ζs)^Q) μt))^(5/2)) (((((ζs)^4)+((ζs)^(2 Q)) (-1+μt)-((ζs)^(3+Q)) μt))^2) ((ρs)^4))};

orcusp2=4;

Assuming[0<ζ<ζs<1&&Q>2&&μt>0&&ρs>0,(equationdimless/.{ρt->Function[ζ,ρs+((∑_(kk=1)^orcusp2)cR[kk](((ζs-ζ))^(kk/2))+((O[ζ,ζs])^((orcusp2+1)/2)))]}/.solcusp2//Simplify)/.((σ2)^(n_)):>If[EvenQ[n],1,σ2]//Simplify]
Solve[%==0,cR[orcusp2]];
Assuming[0<ζ<ζs<1&&Q>2,%//Simplify]

Interestingly, now the sign matters. The two choices will give us two different branches.

Clear[ρtcust2p,ρtcust2m,dρtcust2p,dρtcust2m];

ρtcust2p[Q_,μt_,ρs_,ζs_,ζ_]=ρs+((∑_(kk=1)^orcusp2)cR[kk](((ζs-ζ))^(kk/2)))/.solcusp2/.σ2->1;
dρtcust2p[Q_,μt_,ρs_,ζs_,ζ_]=D[%,ζ];

ρtcust2m[Q_,μt_,ρs_,ζs_,ζ_]=ρs+((∑_(kk=1)^orcusp2)cR[kk](((ζs-ζ))^(kk/2)))/.solcusp2/.σ2->-1;
dρtcust2m[Q_,μt_,ρs_,ζs_,ζ_]=D[%,ζ];

computeSolMutual[Q_,μt_,ρs_,ζs_]:=(solnumericalMutualplus[Q,μt,ρs,ζs]=NDSolve[{eqNum[Q,μt]==0,ρt[ζs-ϵIRnow]==ρtcust2p[Q,μt,ρs,ζs,ζs-ϵIRnow],ρt'[ζs-ϵIRnow]==dρtcust2p[Q,μt,ρs,ζs,ζs-ϵIRnow]},ρt,{ζ,ζs-ϵIRnow,ϵUVnow},WorkingPrecision->30,Method->"StiffnessSwitching"]//Flatten;
solnumericalMutualminus[Q,μt,ρs,ζs]=NDSolve[{eqNum[Q,μt]==0,ρt[ζs-ϵIRnow]==ρtcust2m[Q,μt,ρs,ζs,ζs-ϵIRnow],ρt'[ζs-ϵIRnow]==dρtcust2m[Q,μt,ρs,ζs,ζs-ϵIRnow]},ρt,{ζ,ζs-ϵIRnow,ϵUVnow},WorkingPrecision->30,Method->"StiffnessSwitching"]//Flatten;
)

computeSolMutual[5,0,1/2,1/2]

plotmutual[Q_,μt_,ρs_,ζs_]:=(
If[solnumericalMutualplus[Q,μt,ρs,ζs][[0]]=!= List,Print["solution not found"];Return[];];
Show[plotSol[solnumericalMutualplus[Q,μt,ρs,ζs]],plotSol[solnumericalMutualminus[Q,μt,ρs,ζs]]]
)

computeSolMutual[5,0,2/10,3/4]
plotmutual[5,0,2/10,3/4]

Boundary conditions: Embedding for large spheres

In this case, the functions for the coefficients can be solved even in a more straight forward way, in one step (meaning, we can solve several orders simultaneously).

orIR=5;
equationdimless/.ρt->Function[ζ,(∑_(k=0)^orIR)cR[k](((1-ζ))^k)+((O[ζ,1])^(orIR+1))];
solDiskIR=Solve[%[[3]]==0,Table[cR[iter],{iter,1,orIR}]]//Simplify//Flatten;


Clear[ρtIR,dρtIR]
ρtIR[Q_,μt_,ρIR_,ζ_]=ρt[ζ]/.ρt->Function[ζ,(∑_(k=0)^orIR)cR[k](((1-ζ))^k)]/.solDiskIR/.cR[0]->ρIR//Simplify;
dρtIR[Q_,μt_,ρIR_,ζ_]=D[%,ζ];

Clear[computeSolIR]
computeSolIR[Q_,μt_,ρIR_]:=(solnumericalIR[Q,μt,ρIR]=NDSolve[{eqNum[Q,μt]==0,ρt[1-ϵIRnow]==ρtIR[Q,μt,ρIR,1-ϵIRnow],ρt'[1-ϵIRnow]==dρtIR[Q,μt,ρIR,1-ϵIRnow]},ρt,{ζ,1-ϵIRnow,ϵUVnow},WorkingPrecision->wpc]//Flatten)

Block[{μt=0},computeSolIR[5,0,2/10];]

plotSol[solnumericalIR[5,0,2/10]]

SetDirectory[NotebookDirectory[]];
Export["connected.pdf",plotSol[solnumerical[5,0,1/2]]]
Export["mutual.pdf",plotmutual[5,0,2/10,3/4]]
Export["disconnected.pdf",plotSol[solnumericalIR[5,0,2/10]]]

Computation of entanglement entropies in the particular case Q=5.

To get concrete results it is convenient to fix the parameters. In particular, the independent modes close to the boundary depend on the choice of parameters.

equationdimlessFIX=equationdimless/.Q->5/.σt-> 1/.μt->0//Simplify;

Boundary conditions: Embedding for small spheres

orcusp1=5;
solcusp1={cR[0]->((2 (Sqrt[ζs]))/(Sqrt[3])),cR[1]->((5+4 ((ζs)^6))/(8 (Sqrt[3]) (Sqrt[ζs]) (-1+((ζs)^6)))),cR[2]->((-89+7672 ((ζs)^6)+112 ((ζs)^12))/(2304 (Sqrt[3]) ((ζs)^(3/2)) (((-1+((ζs)^6)))^2))),cR[3]->((-307+167244 ((ζs)^6)+290160 ((ζs)^12)-1472 ((ζs)^18))/(36864 (Sqrt[3]) ((ζs)^(5/2)) (((-1+((ζs)^6)))^3))),cR[4]->(3124211+1884920368 ((ζs)^6)+17413851936 ((ζs)^12)+8019872512 ((ζs)^18)-7204096 ((ζs)^24))/(530841600 (Sqrt[3]) ((ζs)^(7/2)) (((-1+((ζs)^6)))^4)),cR[5]->(-2211607+4647266020 ((ζs)^6)+187777554336 ((ζs)^12)+384507582592 ((ζs)^18)+74771934464 ((ζs)^24)+22139904 ((ζs)^30))/(2831155200 (Sqrt[3]) ((ζs)^(9/2)) (((-1+((ζs)^6)))^5))};

Assuming[0<ζ<ζs<1&&Q>2&&μt>0,(equationdimlessFIX/.{ρt->Function[ζ,(∑_(kk=0)^orcusp1)cR[kk](((ζs-ζ))^(kk+(1/2)))+((O[ζ,ζs])^(orcusp1+(3/2)))]}/.solcusp1//Simplify)/.((σ1)^(n_)):>If[EvenQ[n],1,σ1]//Simplify]
Solve[%==0,cR[orcusp1]];
Assuming[0<ζ<ζs<1,%//Simplify]

Clear[ρtcust1FIX,dρtcust1FIX];
ρtcust1FIX[ζs_,ζ_]=(∑_(kk=0)^orcusp1)cR[kk](((ζs-ζ))^(kk+(1/2)))/.solcusp1/.σ1->1;
dρtcust1FIX[ζs_,ζ_]=D[%,ζ];

computeSolFIX[ζs_]:=(solnumerical[ζs]=NDSolve[{equationdimlessFIX==0,ρt[ζs-ϵIRnow]==ρtcust1FIX[ζs,ζs-ϵIRnow],ρt'[ζs-ϵIRnow]==dρtcust1FIX[ζs,ζs-ϵIRnow]},ρt,{ζ,ζs-ϵIRnow,ϵUVnow},WorkingPrecision->wpc,Method->"StiffnessSwitching"]//Flatten);

computeSolFIX[1/2];

plotSol[solnumerical[1/2]]

Boundary conditions: Embedding for large spheres

orIR=10;
equationdimlessFIX/.ρt->Function[ζ,(∑_(k=0)^orIR)cR[k](((1-ζ))^k)+((O[ζ,1])^(orIR+1))];
solDiskIR=Solve[%[[3]]==0,Table[cR[iter],{iter,1,orIR}]]//Simplify//Flatten;


Clear[ρtIRFIX,dρtIRFIX]
ρtIRFIX[ρIR_,ζ_]=ρt[ζ]/.ρt->Function[ζ,(∑_(k=0)^orIR)cR[k](((1-ζ))^k)]/.solDiskIR/.cR[0]->ρIR//Simplify;
dρtIRFIX[ρIR_,ζ_]=D[%,ζ];

Clear[computeSolIRFIX]
computeSolIRFIX[ρIR_]:=(solnumericalIR[ρIR]=NDSolve[{equationdimlessFIX==0,ρt[1-ϵIRnow]==ρtIRFIX[ρIR,1-ϵIRnow],ρt'[1-ϵIRnow]==dρtIRFIX[ρIR,1-ϵIRnow]},ρt,{ζ,1-ϵIRnow,ϵUVnow},WorkingPrecision->wpc]//Flatten)

Block[{μt=0},computeSolIRFIX[2/10];]

plotSol[solnumericalIR[2/10]]

Boundary conditions: Near the boundary expansion

Near the boundary, the kind of series expansions that we have to consider for the embedding is

ansatzUV[orUV_]:={ρt->Function[ζ,(∑_(k=0)^orUV)(cR[k,0]+cR[k,1]Log[ζ]+cR[k,2]((Log[ζ])^2)+cR[k,3]((Log[ζ])^3))(ζ^k)+((O[ζ,0])^(orUV+1))/.soldiskUV/.cR[0,0]-> a00/.cR[4,0]->a40]}

orUV=1;soldiskUV={};

soldiskUV={cR[0,1]->0,cR[1,0]->0,cR[1,1]->0,cR[2,0]->-(1/(4 cR[0,0])),cR[2,1]->0,cR[3,0]->0,cR[3,1]->0,cR[4,1]->(1/(32 ((cR[0,0])^3))),cR[5,0]->0,cR[5,1]->0,cR[6,0]->((31-384 μt ((cR[0,0])^4)+1536 ((cR[0,0])^3) cR[4,0])/(4608 ((cR[0,0])^5))),cR[6,1]->(1/(96 ((cR[0,0])^5))),
cR[7,0]->0,cR[7,1]->0,cR[7,2]->0,cR[8,0]->(1/(1179648 ((cR[0,0])^7)))(2587+31488 ((cR[0,0])^3) cR[4,0]+768 μt ((cR[0,0])^4) (-11+96 ((cR[0,0])^2)+768 ((cR[0,0])^3) cR[4,0])-73728 ((cR[0,0])^6) (1+28 ((cR[4,0])^2))),cR[8,1]->((41+768 μt ((cR[0,0])^4)-5376 ((cR[0,0])^3) cR[4,0])/(49152 ((cR[0,0])^7))),cR[8,2]->-(7/(4096 ((cR[0,0])^7))),cR[9,0]->0,cR[9,1]->0,cR[9,2]->0,cR[10,0]->(2370041/(8847360000 ((cR[0,0])^9)))+((157 μt)/(30720 ((cR[0,0])^5)))-(((μt)^2)/(20 cR[0,0]))-((73243 cR[4,0])/(2304000 ((cR[0,0])^6)))+((3 μt cR[4,0])/(8 ((cR[0,0])^2)))+((13 (-1+μt-194 ((cR[4,0])^2)))/(2400 ((cR[0,0])^3)))+(2/5) cR[4,0] (1-μt+8 ((cR[4,0])^2)),cR[10,1]->-(1/(73728000 ((cR[0,0])^9)))(73243+57600 μt ((cR[0,0])^4) (-15+16 ((cR[0,0])^2))+4842240 ((cR[0,0])^3) cR[4,0]-921600 ((cR[0,0])^6) (1+24 ((cR[4,0])^2))),cR[10,2]->((-1261+11520 ((cR[0,0])^3) cR[4,0])/(1228800 ((cR[0,0])^9))),cR[10,3]->(1/(10240 ((cR[0,0])^9)))};

soldiskUV=Join[soldiskUV,Table[cR[k,2]->0,{k,0,6}],Table[cR[k,3]->0,{k,0,9}]];

The previous solve the equation up to a sufficiently high order already.

orUV=10;
Block[{Q=5},
equationdimless/.ansatzUV[orUV]]

It is interesting to see what the solution looks like

ρt[ζ]/.ansatzUV[5]

Now the integral that we have to perform is

Clear[ℒ]
 ((dζdR )^(-1)) R ρ[R](Sqrt[1+(R^4)f[R]((ρ'[R])^2)])/.changeRtoζ/.rulef/.R->Rconf /ζ/.ρ->((((Rconf))^(-1))ρt[#]&)/.rulef/.Q->5//Simplify
ℒ[ζ_]=(Rconf^(-1))%

Note the scaling as ℒ ∝(Subscript[R, conf]), that we are factoring out.

ℒ[ζ]/.ansatzUV[6]//Simplify

Then the counterterms that we have to consider are:

ℒct[ζ_]=-(a00/(ζ^3))+(1/(8 a00 ζ));
Sct[a00_,ζ_]=Integrate[ℒct[ζ],ζ]//Simplify

Sct[a00,ζ]

Here we fic μ to zero, otherwhise the integrand would be different

integrand[a00_,ζ_]=ℒ[ζ]-ℒct[ζ]/.μt->0

To minimize the error, we use the expansion to compute the integral from ϵUV to zero analitically.

Clear[integrandUV];
integrandUV[a00_,ζ_]=integrand[a00,ζ]/.ansatzUV[orUV]/.μt->0//Normal;
Integrate[%,ζ];
SUV[a00_,a40_,ζ_]=%-(Limit[%,ζ-> 0])//Simplify

We can compute the subleading coefficients from the near value of the following combinations:

Limit[{ρt[ζ],-((13 ζ+12 ζ Log[ζ]-16 ((ρt[ζ])^3) Derivative[3][ρt][ζ])/(384 ζ ((ρt[ζ])^3)))}/.ansatzUV[4]//Simplify,ζ->0]

Clear[computeEverything,computeEverythingIR]
computeEverything[ζs_,compute_:True]:=Block[
{a00now,a40now,EEnow},

If[compute,computeSolFIX[ζs];];

{a00now,a40now}={ρt[ζ],-((13 ζ+12 ζ Log[ζ]-16 ((ρt[ζ])^3) Derivative[3][ρt][ζ])/(384 ζ ((ρt[ζ])^3)))}/.ζ->ϵUVnow/.solnumerical[ζs];

EEnow=NIntegrate[integrand[a00now,ζ]/.solnumerical[ζs],{ζ,solnumerical[ζs][[1,2,1,1,2]],solnumerical[ζs][[1,2,1,1,1]]},WorkingPrecision->wpc,AccuracyGoal->acc]+ SUV[a00now,a40now,solnumerical[ζs][[1,2,1,1,1]]]-Sct[a00now,ζs];
{a00now,EEnow,ζs,a40now}
]

computeEverythingIR[ρIRnow_,compute_:True]:=Block[
{a00now,a40now,EEnow},

If[compute,computeSolIRFIX[ρIRnow];];

{a00now,a40now}={ρt[ζ],-((13 ζ+12 ζ Log[ζ]-16 ((ρt[ζ])^3) Derivative[3][ρt][ζ])/(384 ζ ((ρt[ζ])^3)))}/.ζ->ϵUVnow/.solnumericalIR[ρIRnow];


EEnow=NIntegrate[integrand[a00now,ζ]/.solnumericalIR[ρIRnow],{ζ,solnumericalIR[ρIRnow][[1,2,1,1,2]],solnumericalIR[ρIRnow][[1,2,1,1,1]]},WorkingPrecision->wpc,AccuracyGoal->acc]-(a00now/2);

{a00now,EEnow,ρIRnow,a40now}

]

Examples:

ϵUVnow=(10^(-5));

computeEverything[1/2,True]
plotSol[solnumerical[1/2]]

computeEverythingIR[1/2,True]
plotSol[solnumericalIR[1/2]]

Full computation

The following computations can be parallelized, if you wish. Otherwise, replace “ParallelTable” by “Table” (it will take quite some time.)

LaunchKernels[8]; (* Choose the number of kernels to do the computations *)

tableEE=ParallelTable[$MinPrecision=wpc; computeEverything[ζs],{ζs,1/100,99/100,1/500}];
tableEE2=ParallelTable[$MinPrecision=wpc;computeEverything[ζs],{ζs,9901/10000,9999/10000,1/50000}];

tableEEIR=
Join[
ParallelTable[$MinPrecision=wpc; computeEverythingIR[ρIR,True],{ρIR,1/1000,1/100,1/2000}],
ParallelTable[$MinPrecision=wpc; computeEverythingIR[ρIR,True],{ρIR,2/100,2,1/200}],
ParallelTable[$MinPrecision=wpc; computeEverythingIR[ρIR,True],{ρIR,21/10,3,1/20}]
];

more=ParallelTable[$MinPrecision=wpc; computeEverythingIR[ρIR,True],{ρIR,31/10,5,1/10}];

tableEEIR=Join[tableEEIR,more]//DeleteDuplicates;

ListLinePlot[
{Join[tableEE,tableEE2][[;;,{1,2}]],tableEEIR[[;;,{1,2}]]},ImageSize->350,
LabelStyle->labelStyle,Frame->True,AspectRatio->(GoldenRatio^(-1)),
PlotStyle->{{Darker[Blue],Thick},{Orange,Thick}},
PlotRangePadding->0,
ImagePadding->{{45,15},{15,5}}
]
ListLinePlot[
{Join[tableEE,tableEE2][[;;,{1,2}]],tableEEIR[[;;,{1,2}]]},ImageSize->350,
LabelStyle->labelStyle,Frame->True,AspectRatio->(GoldenRatio^(-1)),
PlotStyle->{{Darker[Blue],Thick},{Orange,Thick}},
PlotRange->{{.6,.7},{-.2,-.25}},
PlotRangePadding->0,
ImagePadding->{{45,15},{15,5}}
]

To find the point where they creoss we can do:

poscusp1=Do[If[tableEE[[kk,1]]>tableEE[[kk+1,1]],Return[kk];],{kk,1,Length[tableEE]-1}]
poscusp2=Do[If[tableEEIR[[kk,1]]<tableEEIR[[kk+1,1]],Return[kk];],{kk,20,Length[tableEEIR]-1}]
EEofA1=Interpolation[tableEE[[1;;poscusp1,{1,2}]]];
EEofA2=Interpolation[tableEEIR[[poscusp2;;-1,{1,2}]]];
Acritial=A/.FindRoot[EEofA1[A]-EEofA2[A],{A,EEofA1[[1,1,2]],EEofA2[[1,1,1]]}]

c- function (3D)

And now there are two ways to write down the c-functions:

CLM1[A_]= A D[EEofA1[A],A]-EEofA1[A];
CLM2[A_]= A D[EEofA2[A],A]-EEofA2[A];

tableEEall=Join[tableEE,tableEE2][[1;;-1]];
EEofζs1=Interpolation[tableEEall[[;;,{3,2}]],InterpolationOrder->3];
Aofζs1=Interpolation[tableEEall[[;;,{3,1}]],InterpolationOrder->3];

CLM1ofζs[ζs_]= Aofζs1[ζs] ((D[EEofζs1[ζs],ζs])/(D[Aofζs1[ζs],ζs]))-EEofζs1[ζs];


EEofρIR2=Interpolation[tableEEIR[[;;,{3,2}]],InterpolationOrder->3];
AofρIR2=Interpolation[tableEEIR[[;;,{3,1}]],InterpolationOrder->3];
CLM2ofρIR[ρIR_]= AofρIR2[ρIR] ((D[EEofρIR2[ρIR],ρIR])/(D[AofρIR2[ρIR],ρIR]))-EEofρIR2[ρIR];

Show[ParametricPlot[{Aofζs1[ζs],CLM1ofζs[ζs]},{ζs,tableEEall[[1,3]],tableEEall[[-1,3]]},AspectRatio->(GoldenRatio^(-1))
,PlotStyle->{{Darker[Blue],Thickness[.002]}}
],
ParametricPlot[{AofρIR2[ζs],CLM2ofρIR[ζs]},{ζs,tableEEIR[[3,3]],tableEEIR[[-1,3]]},AspectRatio->(GoldenRatio^(-1)),
PlotStyle->{{Orange,Thickness[.002]}}
],
Show[Plot[CLM1[A],{A,EEofA1[[1,1,1]],Acritial},PlotRange->All,
PlotStyle->{{Darker[Blue],Thick}}],
Plot[CLM2[A],{A,Acritial,EEofA2[[1,1,2]]},PlotRange->All,PlotStyle->{{Orange,Thick}}]],

ListLinePlot[{{Acritial,CLM1[Acritial]},{Acritial,CLM2[Acritial]}},PlotStyle->{Black},PlotMarkers->{Automatic, 6}],

PlotRange->{{0,3},{-1,2/10}},
LabelStyle->labelStyle,Frame->True,AspectRatio->1/3,
PlotRangePadding->0,
ImagePadding->{{36,15},{15,5}},
ImageSize->600
]

CLM1ζUV[A_,ϵUV_]=CLM1[A]-((Log[ϵUV])/(4 A));
CLM2ζUV[A_,ϵUV_]=CLM2[A]-((Log[ϵUV])/(4 A));

Block[{ϵUV=1/100},

Show@Table[

Show[


Show[Plot[CLM1ζUV[A,ϵUV]-(As/A)(CLM1ζUV[As,ϵUV]),{A,EEofA1[[1,1,1]],Acritial},PlotRange->All,
PlotStyle->{{Darker[Blue],Thick}}],
Plot[CLM2ζUV[A,ϵUV]-(As/A)(CLM1ζUV[As,ϵUV]),{A,Acritial,EEofA2[[1,1,2]]},PlotRange->All,PlotStyle->{{Orange,Thick}}]],


ListLinePlot[{{Acritial,CLM1ζUV[Acritial,ϵUV]-(As/Acritial)(CLM1ζUV[As,ϵUV])},{Acritial,CLM2ζUV[Acritial,ϵUV]-(As/Acritial)(CLM1ζUV[As,ϵUV])}},PlotStyle->{Black},PlotMarkers->{Automatic, 6}],


PlotRange->{{0,2.2},{-1,4}},
ImageSize->350,
LabelStyle->labelStyle,Frame->True,AspectRatio->(GoldenRatio^(-1)),
PlotStyle->{{Darker[Blue],Thick},{Orange,Thick}},
PlotRangePadding->0,
ImagePadding->{{45,15},{15,5}}


]

,{As,{1/100,1/20,1/4}}]
]

c- function (4D)

Similarly,

EEofA1=Interpolation[tableEE[[1;;poscusp1,{1,2}]]];
EEofA2=Interpolation[tableEEIR[[poscusp2;;-1,{1,2}]]];
Acritial=A/.FindRoot[EEofA1[A]-EEofA2[A],{A,EEofA1[[1,1,2]],EEofA2[[1,1,1]]}]

Ccyl1[A_]= (A D[#,A]+#&)@(A D[EEofA1[A],A]-EEofA1[A]);
Ccyl2[A_]= (A D[#,A]+#&)@(A D[EEofA2[A],A]-EEofA2[A]);

Show[Plot[Ccyl1[A],{A,0.01,Acritial},PlotRange->All,
PlotStyle->{{Darker[Blue],Thick}}],
Plot[Ccyl2[A],{A,Acritial,EEofA2[[1,1,2]]},PlotStyle->{{Orange,Thick}},PlotRange->All],PlotRange->{{0,3},{-1,3}},AxesOrigin->{0,0},


LabelStyle->labelStyle,Frame->True,AspectRatio->1/2,
PlotRangePadding->0,
ImagePadding->{{36,15},{15,5}},
ImageSize->400,AxesOrigin->{0,0}
]

Export["cfunction_cyl.pdf",%]

intorder=6;

EEofρIR2=Interpolation[tableEEIR[[;;,{3,2}]],InterpolationOrder->intorder];
AofρIR2=Interpolation[tableEEIR[[;;,{3,1}]],InterpolationOrder->intorder];


auxfunc2[ρIR_]= AofρIR2[ρIR] ((D[EEofρIR2[ρIR],ρIR])/(D[AofρIR2[ρIR],ρIR]))-EEofρIR2[ρIR];

CLM2cylofζs[ρIR_]=(3/(2π)) (1/2)AofρIR2[ρIR](AofρIR2[ρIR] ((D[auxfunc2[ρIR],ρIR])/(D[AofρIR2[ρIR],ρIR]))+auxfunc2[ρIR]);

tableEEall=Join[tableEE,tableEE2][[1;;-1]];
EEofζs1=Interpolation[tableEEall[[;;,{3,2}]],InterpolationOrder->intorder];
Aofζs1=Interpolation[tableEEall[[;;,{3,1}]],InterpolationOrder->intorder];

auxfunc1[ζs_]= Aofζs1[ζs] ((D[EEofζs1[ζs],ζs])/(D[Aofζs1[ζs],ζs]))-EEofζs1[ζs];

CLM1cylofζs[ζs_]=(3/(2π))(1/2)Aofζs1[ζs]( Aofζs1[ζs] ((D[auxfunc1[ζs],ζs])/(D[Aofζs1[ζs],ζs]))+auxfunc1[ζs]);
tableEEall=Join[tableEE,tableEE2][[1;;-1]];
EEofζs1=Interpolation[tableEEall[[;;,{3,2}]],InterpolationOrder->intorder];
Aofζs1=Interpolation[tableEEall[[;;,{3,1}]],InterpolationOrder->intorder];

auxfunc1[ζs_]= Aofζs1[ζs] ((D[EEofζs1[ζs],ζs])/(D[Aofζs1[ζs],ζs]))-EEofζs1[ζs];

CLM1cylofζs[ζs_]=(3/(2π))(1/2)Aofζs1[ζs]( Aofζs1[ζs] ((D[auxfunc1[ζs],ζs])/(D[Aofζs1[ζs],ζs]))+auxfunc1[ζs]);

ζsturn=ζs/.FindRoot[Aofζs1[ζs]-Acritial,{ζs,tableEEall[[-510,3]]}]
ρIturn=ζs/.FindRoot[AofρIR2[ζs]-Acritial,{ζs,tableEEall[[-510,3]]}]

labelStyle=Directive[,FontFamily->"Courier",FontSize->14];

plot=Show[
ParametricPlot[{Aofζs1[ζs],CLM1cylofζs[ζs]},{ζs,tableEEall[[4,3]],ζsturn},AspectRatio->(GoldenRatio^(-1))
,PlotStyle->{{Darker[Blue],Thick}},PlotRange->All
],

ParametricPlot[{Aofζs1[ζs],CLM1cylofζs[ζs]},{ζs,ζsturn,tableEEall[[-510,3]]},AspectRatio->(GoldenRatio^(-1))
,PlotStyle->{{Darker[Blue],Thickness[.002]}},PlotRange->All
],

ParametricPlot[{AofρIR2[ζs],CLM2cylofζs[ζs]},{ζs,tableEEIR[[55,3]],0.2977679106110384`},AspectRatio->(GoldenRatio^(-1)),
PlotStyle->{{Orange,Thickness[0.002]}},PlotRange->All
],
ParametricPlot[{AofρIR2[ζs],CLM2cylofζs[ζs]},{ζs,0.2977679106110384`,tableEEIR[[-1,3]]},AspectRatio->(GoldenRatio^(-1)),
PlotStyle->{{Orange,Thick}},PlotRange->All
],
ListLinePlot[{{Acritial,CLM2cylofζs[ρIturn]},{Acritial,CLM1cylofζs[ζsturn]}},PlotStyle->Black,PlotMarkers->{Automatic, 5}]
,
PlotRange->{{0,2.5},{-.25,.1}},
ImageSize->350,
LabelStyle->labelStyle,Frame->True,AspectRatio->(GoldenRatio^(-1)),
PlotStyle->{{Darker[Blue],Thick},{Orange,Thick}},
PlotRangePadding->0,
ImagePadding->{{45,15},{15,5}}

(*
PlotRange->{{0,3},{-.5,.3}},
LabelStyle->labelStyle,Frame->True,AspectRatio->1/2,
PlotRangePadding->0,
ImagePadding->{{36,15},{15,5}},
ImageSize->400,AxesOrigin->{0,0}*)
]
Export["cfunction_cyl_new.pdf",%]
Show[plot,PlotRange-> {{.645,.675},{-.22,-0.08}},LabelStyle->Directive[,FontFamily->"Courier",FontSize->12],ImageSize-> 350/2,FrameTicks-> {{{-.1,-.2},None},{{.65,.67},None}}]
Export["inset.pdf",%]