Figure 1

Figure 1a

Quit

(*mass term from equation (1.3) in the letter*)
al=0.1
Bz:=0.05
Plot[((-4*al*(Bz^2))/(ρ^4)),{ρ,0.01,2},PlotRange->{0,-5},PlotStyle->{Gray},PlotLegends->{"k=0"}]

(*The eoms in this section are the Euler Lagrangian equation of the action in equation 1.2 in the letter*)
(*If L=0, eom for a[rho] reduce to:*)
DSolve[a''[ρ]+3((a'[ρ])/ρ)==(k^2)((a[ρ])/(ρ^4)),a[ρ],ρ]
(*solutions are Bessel functions:*)

(*The large ρ asymptotic expansion of the general solution is given as*)
Series[((BesselK[1,((Sqrt[(k^2)])/ρ)] (C[1]))/(2 (Sqrt[(k^2)]) (Sqrt[π]) ρ))+((2 BesselI[1,((Sqrt[(k^2)])/ρ)] (C[2]))/((Sqrt[(k^2)]) ρ)),{ρ,Infinity,2}]

(*To make the aysmptotic expansion at UV satifies B=k*a[k], C2=0, and C1 ∝ k*)

(*The final solution with L=0 has the following form:*)
α:=0.1
B:=0.05
(*k:=0.173*)
afun[ρ_]:=B*((BesselK[1,(k/ρ)])/ρ)
ap[ρ_]:=B*((k BesselK[0,(k/ρ)]-2 ρ BesselK[1,(k/ρ)]+k BesselK[2,(k/ρ)])/(2 (ρ^3)))

k0:=0.173(*parameter k to adjust with k=0.173, 0.218, 0.2, 0.24*)
color={{Blue},{Black,Dashed},Orange,Red}
ps:=color[[4]]
lg:={"k="<>ToString[k]}
Module[{},
k:=k0;
Plot[k*afun[ρ],{ρ,0,100},PlotRange->Full,AxesLabel->{"ρ","k*a(ρ)"},PlotLegends->lg,PlotStyle->ps]]
Module[{},
k:=k0;
Plot[-4α((k^2)(((afun[ρ])^2)/(ρ^4))+((ap[ρ])^2)),{ρ,0.0,5},PlotRange->{{0,1},{-2.6,0}},AxesLabel->{"ρ","(Subscript[m, L]^2)"},PlotLegends->lg,PlotStyle->ps]]

Show[%499,%507,%515,%525,PlotRange->{{0,1},{0,0.05}}]
Show[%500,%508,%516,%526,PlotRange->{{0,1},{-2.5,0}},AxesLabel->{"ρ","(Subscript[m, L]^2)"}]

Show[%544,%549](*figure 1*)

Figure 1b

ClearAll["Global`*"];

(*X={t,x,y,z,ρ}*)
ξ={t,x,y,z,ρ};
Gdd=DiagonalMatrix@{-(ρ^2),(ρ^2),(ρ^2),(ρ^2),(1/(ρ^2))}
Avec={0,a[ρ] Cos[k*z],a[ρ] Sin[k*z],0,0};
Fmat=Table[((∂_(ξ[[μ]]))Avec[[ν]]-(∂_(ξ[[ν]]))Avec[[μ]]),{μ,5},{ν,5}]
%//MatrixForm
(*F2 is in fact the (F^(μν))(Subscript[F, μν])*)
GUU=Inverse[Gdd]
FUU=GUU.Fmat.GUU;
%//MatrixForm
F2=Tr[-GUU.Fmat.GUU.Fmat]//FullSimplify

(*Plug the flutuation expansion into the action*)lag=(1/ρ)(-(1/2)(-((((∂_t)(L0[ρ]+ϵ*L1[t,x,y,z,ρ])))^2)+((((∂_x)(L0[ρ]+ϵ*L1[t,x,y,z,ρ])))^2)+((((∂_y)(L0[ρ]+ϵ*L1[t,x,y,z,ρ])))^2)+((((∂_z)(L0[ρ]+ϵ*L1[t,x,y,z,ρ])))^2)))+(ρ^3)(-(1/2)((((∂_ρ)(L0[ρ]+ϵ*L1[t,x,y,z,ρ])))^2))-(1/4)(ρ^3)F2+α ρ F2 (((L0[ρ]+ϵ*L1[t,x,y,z,ρ]))^2)
(*effective action for the fluctuation*)
lagfluc=SeriesCoefficient[%,{ϵ,0,2}]

Needs["VariationalMethods`"]
EulerEquations[lagfluc,{ L1[t,x,y,z,ρ]},{t,x,y,z,ρ}]

(*Plane wave expansion L1=L0[ρ]*Exp[ik.x], gave the following equation of motion*)
(ρ^6) L0''[ρ]+3 (ρ^5) L0'[ρ]+(4*(k^2)*alpha*((a[ρ])^2)+4*alpha*(ρ^4)*((aprime[ρ])^2)+M2*(ρ^2))*L0[ρ]==0;

solve the fluctuation of L numerically

plot the pole to search the mass of the fluctuations

(*---1. Define Parameters and Functions---*)
alpha=0.1;
k=0.1712;
ir=0.01; (*Initial rho value*)
fr=100;  (*Final rho value*)

(*Define the function a[rho] and its derivative a'[rho]*)
a[ρ_]:=(0.05*((BesselK[1,((Abs[k])/ρ)])/ρ));
(*aprime[ρ]=D[a[ρ]]*)
aprime[ρ_]:=(1/(ρ^3))(0.025` Abs[k] BesselK[0,((Abs[k])/ρ)]-0.05` ρ BesselK[1,((Abs[k])/ρ)]+0.025` Abs[k] BesselK[2,((Abs[k])/ρ)]);

(*---2. Create a Function to Solve the ODE for a given M2---*)
(*This function returns the value of the solution at the endpoint fr*)
endValue[m2_?NumericQ]:=Module[{ode,ic},
M2=m2;
ode=(ρ^6) L0''[ρ]+3 (ρ^5) L0'[ρ]+(4*(k^2)*alpha*((a[ρ])^2)+4*alpha*(ρ^4)*((aprime[ρ])^2)+M2*(ρ^2))*L0[ρ]==0;(*the equation of motion for the fluctuation L0*)
ic={L0[ir]==1,L0'[ir]==0};
NDSolveValue[{ode,ic},L0[fr],{ρ,ir,fr}]]

(*---3. Scan M2 and Generate Data---*)
scanData=Table[{m2,endValue[m2]},{m2,-0.001,0.006,0.00002}];


(*---4. Plot the Scan Results---*)
ListLinePlot[scanData,PlotRange->All,AxesLabel->{"M2","1/|L0[fr]|"},PlotLabel->"Scan of M2 vs. L0[fr]",GridLines->Automatic,Epilog->{Red,Dashed,(*Add a vertical line at the value of M2 found previously*)Line[{{0.00147656,0},{0.00147656,Max[scanData[[All,2]]]}}]}]

fine-tune to find the mass square of the meson

(*For each k and masse square saved in the data section below, can be setted up into the parameter k and m2 guess to get the output of the plot*)

ClearAll["Global`*"]
(*---1. Define Parameters and Functions---*)
alpha=0.1;
k=0.1712;(*set up the k value*)
m2guess=0.005(*make a guess of the mass square *)
ir=0.01; (*IR rho value*)
fr=100;  (*UV rho value*)


(*Define the function a[rho] and its derivative a'[rho]*)
a[ρ_]:=(0.05*((BesselK[1,((Abs[k])/ρ)])/ρ));
(*aprime[ρ]=D[a[ρ]]*)
aprime[ρ_]:=(1/(ρ^3))(0.025` Abs[k] BesselK[0,((Abs[k])/ρ)]-0.05` ρ BesselK[1,((Abs[k])/ρ)]+0.025` Abs[k] BesselK[2,((Abs[k])/ρ)]);

(*---2. Create a Function for the Root Finder---*)
(*This function takes a numerical value for M2,solves the ODE,and returns the value at the endpoint L0[fr]*)
endValue[m2_?NumericQ]:=Module[{sol,ode,ic},
M2=m2;
ode=(ρ^6) L0''[ρ]+3 (ρ^5) L0'[ρ]+(4*(k^2)*alpha*((a[ρ])^2)+4*alpha*(ρ^4)*((aprime[ρ])^2)+M2*(ρ^2))*L0[ρ]==0;
ic={L0[ir]==1,L0'[ir]==0};
sol=NDSolveValue[{ode,ic},L0[fr],{ρ,ir,fr}];
sol]

(*---3. Find the Root---*)
(*Use FindRoot to find the value of m2 for which endValue[m2] is 0. Start with an initial guess of 0.0001.*)
m2Solution=FindRoot[endValue[m2],{m2,m2guess}];
M2=m2/. m2Solution;

(*Print the found value of M2*)
Print["The required value for M2 is: ",M2];


(*---4. Solve the ODE Numerically with the new M2---*)
odeFinal=(ρ^6) L0''[ρ]+3 (ρ^5) L0'[ρ]+(4*(k^2)*alpha*((a[ρ])^2)+4*alpha*(ρ^4)*((aprime[ρ])^2)+M2*(ρ^2))*L0[ρ]==0;
icFinal={L0[ir]==1,L0'[ir]==0};
solFinal=NDSolve[{odeFinal,icFinal},L0,{ρ,ir,fr}];


(*---5. Plot the Result---*)
(*Visualize the final solution L0[rho] as a function of rho*)
Plot[Evaluate[L0[ρ]/. solFinal],{ρ,ir,100},PlotRange->All,AxesLabel->{"ρ","L0[ρ]"},PlotLabel->"Solution of the ODE with M2 = "<>ToString[M2],Epilog->{Red,PointSize[Large],Point[{fr,0}]}]
(*---5. Plot the Result---*)
(*Visualize the final solution L0[rho] as a function of rho*)
Plot[Evaluate[L0[ρ]/. solFinal],{ρ,ir,2},PlotRange->All,AxesLabel->{"ρ","L0[ρ]"},PlotLabel->"Solution of the ODE with M2 = "<>ToString[M2],Epilog->{Red,PointSize[Large],Point[{fr,0}]}]

data

(*data for k=0.17*)
k=0.17
(*excitation number*)
n={0,1,2}
(*mass square*)
m2={-0.00009400725781568038`,0.0011341552816260495`,0.004907053139649996`}

k=0.1712
n={0,1,2}
m2={3.1402898254930186`*^-6,0.001197663985276833`,0.0049626546214579514`}

k=0.173;
n=  {0,1,2}
m2={0.00012077857060824848`,0.001302331603818405`,0.005047801345771665`}

k=0.1733;
n=  {0,1,2}
m2={0.000137457,0.0013206122600509496`,0.005062170595062094`}

(*re-organise the data*)
(*{{k, meson mass}} with varied excited number n*)
(*n=0*)
datan0={{0.17,-0.00009400725781568038},{0.1712,3.1402898254930186*^-6},{0.173,0.00012077857060824848},{0.1733,0.000137457}};

(*n=1*)
datan1={{0.17,0.0011341552816260495},{0.1712,0.001197663985276833},{0.173,0.001302331603818405},{0.1733,0.0013206122600509496}};

(*n=2*)
datan2={{0.17,0.004907053139649996},{0.1712,0.0049626546214579514},{0.173,0.005047801345771665},{0.1733,0.005062170595062094}};

ListLinePlot[{datan0,datan1,datan2},AxesLabel->{"k",None},PlotStyle->{Red,Green,Blue},PlotLegends->{"n = 0","n = 1","n = 2"}]

Figure 2

k=0 (Bz=0.05)

at k= 0, the EoM eq.(2.4) reduce to

ClearAll["Global`*"]
{L'[ρ] (-4 (Bz^2) β-3 (ρ^4)-4 (k^2) β ((a[ρ])^2)+8 (k^2) β ρ a[ρ] a'[ρ]+4 β (ρ^4) a'[ρ] (3 a'[ρ]+2 ρ a''[ρ]))+ρ (4 (Bz^2) β-(ρ^4)+4 (k^2) β ((a[ρ])^2)+4 β (ρ^4) ((a'[ρ])^2)) L''[ρ]==((2 L[ρ] (2 (Bz^2) α-γ (ρ^2)+2 (k^2) α ((a[ρ])^2)+2 α (ρ^4) ((a'[ρ])^2)))/ρ),(((4 α ((L[ρ])^2)+(ρ^2) (-1+4 β ((L'[ρ])^2))) ((k^2) a[ρ]-(ρ^4) a''[ρ]))/(ρ^3))==a'[ρ] (4 α ((L[ρ])^2)+8 α ρ L[ρ] L'[ρ]+(ρ^2) (-3+4 β L'[ρ] (3 L'[ρ]+2 ρ L''[ρ])))}/.{(*Bz->0,*)β->0,γ->0,k->0}
Solve[%,{L''[ρ],a''[ρ]}][[1]]

(*in the limit of k=0, a[ρ] approaches to a constant solution, and the derivative a'[ρ] approaches to zero *)
L''[ρ]->((-4 (Bz^2) α L[ρ]-4 α (ρ^4) L[ρ] ((a'[ρ])^2)-3 (ρ^5) L'[ρ])/(ρ^6))/.{a'[ρ]->0}

(*---Part 1:Re-solve the system with the found parameters---*)
ClearAll["Global`*"];
finalrir=0.10032867197814097`;
α:=0.1
Bz:=0.05
uv:=100
finalSolution=NDSolve[{-3 (ρ^4) L'[ρ]-(ρ^5) L''[ρ]==((2 L[ρ] (2 (Bz^2) α))/ρ),L[finalrir]==finalrir,L'[finalrir]==0},{L,a},{ρ,finalrir,uv}];

(*Get the final interpolating functions*)
LsolFinal=L/. finalSolution[[1]];
(*asolFinal=a/. finalSolution[[1]];*)

(*Plot L(ρ) and k*a(ρ) on the same graph for comparison*)
plot1=Plot[LsolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","L[ρ]"},PlotRange->{{0,10},Full},PlotStyle->Gray,PlotLegends->{"k=0"}]
(*plota=Plot[asolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","k a[ρ]"}]*)

k=0.05

(*Clear definitions*)Clear[α,k,B,uv,SolveWithParameter]

(*Define constants with higher precision by using backticks ``*)
α=1/10;
k=5/100;
g1=0.1;
g2=0.047/0.05;
uv=100;
B=5/100;

(*Set the desired working precision*)
wp=30;

(*Function to solve the ODEs and return the residuals*)
SolveWithParameter[ir_?NumericQ,a0_?NumericQ]:=Module[{solution,Lsol,asol,L,a,ρ},solution=Quiet@NDSolve[{L''[ρ]==(-2 L[ρ] (2 k^2 α a[ρ]^2+2 α ρ^4 a'[ρ]^2)-3 ρ^5 L'[ρ])/ρ^6,a''[ρ]==-(((-ρ^4 (k^2 a[ρ] (ρ^2-4 α L[ρ]^2)+ρ^3 (-3 ρ^2+4 α L[ρ]^2) a'[ρ])-8 α ρ^8 L[ρ] a'[ρ] L'[ρ])/(ρ^8 (ρ^2-4 α L[ρ]^2)))),L[ir]==ir,L'[ir]==0,a[ir]==a0,a'[ir]==0},{L,a},{ρ,ir,uv},Method->"StiffnessSwitching",WorkingPrecision->wp,AccuracyGoal->wp-10,PrecisionGoal->wp-10];
Lsol=L/. solution[[1]];
asol=a/. solution[[1]];
{asol[uv]*k-B,Lsol[uv]}]

(*Find the root using the corrected function name and higher precision*)
foundParameters=Quiet@FindRoot[SolveWithParameter[rir,a0],{{rir,g1},{a0,g2}},MaxIterations->2000,WorkingPrecision->wp];

(*Display the found parameters*)
Print["Found parameters:"];
foundParameters

color=RGBColor[0.6,0.8,0.4]
legend={"k=0.05"}
(*Extract the final values for rir and a0 from the FindRoot result*)
finalrir=rir/. foundParameters
finala0=a0/. foundParameters

(*Solve the ODEs one last time with the optimal parameters to get the solution functions*)
finalSolution=NDSolve[{L''[ρ]==(-2 L[ρ] (2 k^2 α a[ρ]^2+2 α ρ^4 a'[ρ]^2)-3 ρ^5 L'[ρ])/ρ^6,a''[ρ]==-((-ρ^4 (k^2 a[ρ] (ρ^2-4 α L[ρ]^2)+ρ^3 (-3 ρ^2+4 α L[ρ]^2) a'[ρ])-8 α ρ^8 L[ρ] a'[ρ] L'[ρ])/(ρ^8 (ρ^2-4 α L[ρ]^2))),L[finalrir]==finalrir,L'[finalrir]==0,a[finalrir]==finala0,a'[finalrir]==0},{L,a},{ρ,finalrir,uv},Method->"StiffnessSwitching",AccuracyGoal->20,PrecisionGoal->20,WorkingPrecision->30];

(*Get the final interpolating functions*)
LsolFinal=L/. finalSolution[[1]];
asolFinal=a/. finalSolution[[1]];

(*Plot L(ρ) and k*a(ρ) on the same graph for comparison*)
plot1=Plot[LsolFinal[ρ],{ρ,finalrir,uv},PlotStyle->color,PlotLegends->legend,AxesLabel->{"ρ","L[ρ]"},PlotRange->{{0,10},Full}]
plota=Plot[k*asolFinal[ρ],{ρ,finalrir,uv},PlotStyle->color,PlotLegends->legend,AxesLabel->{"ρ","k a[ρ]"},PlotRange->{0,0.06}]
(*Plot the expression k^2*a(ρ)^2/ρ^3*)
plot2=Plot[(-k^2*asolFinal[ρ]^2)/ρ^3,{ρ,finalrir,uv},PlotStyle->color,AxesLabel->{"ρ","-(((k^2) ((a[ρ])^2))/(ρ^3))"},PlotLegends->legend,PlotRange->{-4,0}]

k=0.1

(*Clear definitions*)
Clear[α,k,B,uv,SolveWithParameter]
(*we use corrected eoms, sgn=+1*)
(*Define constants*)
α:=0.1
k:=0.1
g1:=0.089672
g2:=0.38458
uv:=100
B:=0.05
(*Function to solve the ODEs and return the residuals*)

legend={"k="<>ToString[k]}
(**)
SolveWithParameter[ir_?NumericQ,a0_?NumericQ]:=Module[{solution,Lsol,asol},solution=Quiet@NDSolve[{L''[ρ]==(-2 L[ρ] (2 k^2 α a[ρ]^2+2 α ρ^4 a'[ρ]^2)-3 ρ^5 L'[ρ])/ρ^6,a''[ρ]==-((-ρ^4 (k^2 a[ρ] (ρ^2-4 α L[ρ]^2)+ρ^3 (-3 ρ^2+4 α L[ρ]^2) a'[ρ])-8 α ρ^8 L[ρ] a'[ρ] L'[ρ])/(ρ^8 (ρ^2-4 α L[ρ]^2))),L[ir]==ir,L'[ir]==0,a[ir]==a0,a'[ir]==0},{L,a},{ρ,ir,uv},Method->"StiffnessSwitching",WorkingPrecision->35];
Lsol=L/. solution[[1]];
asol=a/. solution[[1]];
{asol[uv]*k-B,Lsol[uv]}]
(*Find the root using the corrected function name*)
foundParameters=FindRoot[SolveWithParameter[rir,a0],{{rir,g1},{a0,g2}},(*MaxIterations->100*)WorkingPrecision->35];
(*Display the found parameters*)
Print["Found parameters:"];
foundParameters

color=RGBColor[0.2,0.6,0.8]
(*Extract the final values for rir and a0 from the FindRoot result*)
finalrir=rir/. foundParameters;
finala0=a0/. foundParameters;

(*Solve the ODEs one last time with the optimal parameters to get the solution functions*)
finalSolution=NDSolve[{L''[ρ]==(-2 L[ρ] (2 k^2 α a[ρ]^2+2 α ρ^4 a'[ρ]^2)-3 ρ^5 L'[ρ])/ρ^6,a''[ρ]==-((-ρ^4 (k^2 a[ρ] (ρ^2-4 α L[ρ]^2)+ρ^3 (-3 ρ^2+4 α L[ρ]^2) a'[ρ])-8 α ρ^8 L[ρ] a'[ρ] L'[ρ])/(ρ^8 (ρ^2-4 α L[ρ]^2))),L[finalrir]==finalrir,L'[finalrir]==0,a[finalrir]==finala0,a'[finalrir]==0},{L,a},{ρ,finalrir,uv}];

(*Get the final interpolating functions*)
LsolFinal=L/. finalSolution[[1]];
asolFinal=a/. finalSolution[[1]];

(*Plot L(ρ) and k*a(ρ) on the same graph for comparison*)
plot1=Plot[LsolFinal[ρ],{ρ,finalrir,uv},PlotStyle->color,PlotLegends->legend,AxesLabel->{"ρ","L[ρ]"},PlotRange->{{0,10},Full}]
plota=Plot[k*asolFinal[ρ],{ρ,finalrir,uv},PlotStyle->color,PlotLegends->legend,AxesLabel->{"ρ","k a[ρ]"},PlotRange->{0,0.06}]
(*Plot the expression k^2*a(ρ)^2/ρ^3*)
plot2=Plot[(-k^2*asolFinal[ρ]^2)/ρ^3,{ρ,finalrir,uv},PlotStyle->color,AxesLabel->{"ρ","-(((k^2) ((a[ρ])^2))/(ρ^3))"},PlotLegends->legend,PlotRange->{-4,0}]

k = 0.1733

solution 1

(*---Part 1:Re-solve the system with the found parameters---*)
ClearAll["Global`*"];
(*Extract the final values for rir and a0 from the FindRoot result*)
finalrir=0.03346;
finala0=0.010593386640674705695;
α:=0.1
k:=0.1733
uv:=100
color={Purple,Dashed};
lg={"k="<>ToString[k]}
(*B:=0.06*)
(*Solve the ODEs one last time with the optimal parameters to get the solution functions*)
finalSolution=NDSolve[{L''[ρ]==(-2 L[ρ] (2 k^2 α a[ρ]^2+2 α ρ^4 a'[ρ]^2)-3 ρ^5 L'[ρ])/ρ^6,a''[ρ]==-((-ρ^4 (k^2 a[ρ] (ρ^2-4 α L[ρ]^2)+ρ^3 (-3 ρ^2+4 α L[ρ]^2) a'[ρ])-8 α ρ^8 L[ρ] a'[ρ] L'[ρ])/(ρ^8 (ρ^2-4 α L[ρ]^2))),L[finalrir]==finalrir,L'[finalrir]==0,a[finalrir]==finala0,a'[finalrir]==0},{L,a},{ρ,finalrir,uv},Method->"StiffnessSwitching"];

(*Get the final interpolating functions*)
LsolFinal=L/. finalSolution[[1]];
asolFinal=a/. finalSolution[[1]];

(*Plot L(ρ) and k*a(ρ) on the same graph for comparison*)
plotL1=Plot[LsolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","L[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color,PlotLegends->lg]
plota1=Plot[k*asolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","k a[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color,PlotLegends->lg]
(*Plot the expression k^2*a(ρ)^2/ρ^3*)
plot2=Plot[-4α((((k^2)*((asolFinal[ρ])^2))/(ρ^4))+((asolFinal'[ρ])^2)),{ρ,finalrir,1},PlotStyle->color,PlotLegends->lg,AxesLabel->{"ρ","-4α((((k^2)*((a[ρ])^2))/(ρ^4))+a'[ρ]^2)"},PlotRange->{{0,1},Full}]

solution  2

{{k=0.1733},{rir->0.041896597739897715520138876462434347,a0->0.026088595172430099595592348009630510}}

(*(*---Part 1:Re-solve the system with the found parameters---*)
ClearAll["Global`*"];*)

color=Green
(*Extract the final values for rir and a0 from the FindRoot result*)
finalrir=0.041896597739897715520138876462434347;
finala0=0.026088595172430099595592348009630510;
α:=0.1
k:=0.1733
uv:=100
color=Purple;
(*B:=0.06*)
(*Solve the ODEs one last time with the optimal parameters to get the solution functions*)
finalSolution=NDSolve[{L''[ρ]==(-2 L[ρ] (2 k^2 α a[ρ]^2+2 α ρ^4 a'[ρ]^2)-3 ρ^5 L'[ρ])/ρ^6,a''[ρ]==-((-ρ^4 (k^2 a[ρ] (ρ^2-4 α L[ρ]^2)+ρ^3 (-3 ρ^2+4 α L[ρ]^2) a'[ρ])-8 α ρ^8 L[ρ] a'[ρ] L'[ρ])/(ρ^8 (ρ^2-4 α L[ρ]^2))),L[finalrir]==finalrir,L'[finalrir]==0,a[finalrir]==finala0,a'[finalrir]==0},{L,a},{ρ,finalrir,uv},Method->"StiffnessSwitching"];

(*Get the final interpolating functions*)
LsolFinal=L/. finalSolution[[1]];
asolFinal=a/. finalSolution[[1]];

(*Plot L(ρ) and k*a(ρ) on the same graph for comparison*)
plotL2=Plot[LsolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","L[ρ]"},PlotRange->{{0,2},Full},PlotStyle->color]
plota2=Plot[k*asolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","k a[ρ]"},PlotRange->{{0,2},Full},PlotStyle->color]
(*Plot the expression k^2*a(ρ)^2/ρ^3*)
plot2=Plot[-4α((((k^2)*((asolFinal[ρ])^2))/(ρ^4))+((asolFinal'[ρ])^2)),{ρ,finalrir,1},PlotStyle->color,PlotLegends->{"k=0.1733"},AxesLabel->{"ρ","-4α((((k^2)*((a[ρ])^2))/(ρ^4))+a'[ρ]^2)"},PlotRange->{{0,1},Full}]

finalrir-LsolFinal[finalrir]

Show[Plot[x,{x,0,0.06},PlotStyle->Black],plotL1,plotL2,PlotRange->{{0,1},{0,0.043}},ImageSize->400]
Show[plota1,plota2,PlotRange->{{0,2},{0,0.05}},ImageSize->400]

k = 0.173

solution 1

(*---Part 1:Re-solve the system with the found parameters---*)
ClearAll["Global`*"];
(*Extract the final values for rir and a0 from the FindRoot result*)
finalrir=0.03013196406430698`;
finala0=0.006440761542702392`;
α:=0.1
k:=0.173
uv:=100
color={Red,Dashed};
lg={"k="<>ToString[k]}
(*B:=0.06*)
(*Solve the ODEs one last time with the optimal parameters to get the solution functions*)
finalSolution=NDSolve[{L''[ρ]==(-2 L[ρ] (2 k^2 α a[ρ]^2+2 α ρ^4 a'[ρ]^2)-3 ρ^5 L'[ρ])/ρ^6,a''[ρ]==-((-ρ^4 (k^2 a[ρ] (ρ^2-4 α L[ρ]^2)+ρ^3 (-3 ρ^2+4 α L[ρ]^2) a'[ρ])-8 α ρ^8 L[ρ] a'[ρ] L'[ρ])/(ρ^8 (ρ^2-4 α L[ρ]^2))),L[finalrir]==finalrir,L'[finalrir]==0,a[finalrir]==finala0,a'[finalrir]==0},{L,a},{ρ,finalrir,uv},Method->"StiffnessSwitching"];

(*Get the final interpolating functions*)
LsolFinal=L/. finalSolution[[1]];
asolFinal=a/. finalSolution[[1]];

(*Plot L(ρ) and k*a(ρ) on the same graph for comparison*)
plotL1=Plot[LsolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","L[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color,PlotLegends->lg]
plota1=Plot[k*asolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","k a[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color,PlotLegends->lg]
(*Plot the expression k^2*a(ρ)^2/ρ^3*)
plot2=Plot[-4α((((k^2)*((asolFinal[ρ])^2))/(ρ^4))+((asolFinal'[ρ])^2)),{ρ,finalrir,1},PlotStyle->color,PlotLegends->lg,AxesLabel->{"ρ","-4α((((k^2)*((a[ρ])^2))/(ρ^4))+a'[ρ]^2)"},PlotRange->{{0,1},Full}]

solution  2

(*(*---Part 1:Re-solve the system with the found parameters---*)
ClearAll["Global`*"];*)
(*Extract the final values for rir and a0 from the FindRoot result*)
finalrir=0.04404588987996323`;
finala0=0.03113714045007225`;
α:=0.1
k:=0.173
uv:=100
color={Red};
(*B:=0.06*)
(*Solve the ODEs one last time with the optimal parameters to get the solution functions*)
finalSolution=NDSolve[{L''[ρ]==(-2 L[ρ] (2 k^2 α a[ρ]^2+2 α ρ^4 a'[ρ]^2)-3 ρ^5 L'[ρ])/ρ^6,a''[ρ]==-((-ρ^4 (k^2 a[ρ] (ρ^2-4 α L[ρ]^2)+ρ^3 (-3 ρ^2+4 α L[ρ]^2) a'[ρ])-8 α ρ^8 L[ρ] a'[ρ] L'[ρ])/(ρ^8 (ρ^2-4 α L[ρ]^2))),L[finalrir]==finalrir,L'[finalrir]==0,a[finalrir]==finala0,a'[finalrir]==0},{L,a},{ρ,finalrir,uv},Method->"StiffnessSwitching"];

(*Get the final interpolating functions*)
LsolFinal=L/. finalSolution[[1]];
asolFinal=a/. finalSolution[[1]];

(*Plot L(ρ) and k*a(ρ) on the same graph for comparison*)
plotL2=Plot[LsolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","L[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color]
plota2=Plot[k*asolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","k a[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color]
(*Plot the expression k^2*a(ρ)^2/ρ^3*)
plot2=Plot[-4α((((k^2)*((asolFinal[ρ])^2))/(ρ^4))+((asolFinal'[ρ])^2)),{ρ,finalrir,1},PlotStyle->color,PlotLegends->{"k=0.173"},AxesLabel->{"ρ","-4α((((k^2)*((a[ρ])^2))/(ρ^4))+a'[ρ]^2)"},PlotRange->{{0,1},Full}]

Show[Plot[x,{x,0,0.06},PlotStyle->Black],plotL1,plotL2,PlotRange->{{0,1},{0,0.043}},ImageSize->400]
Show[plota1,plota2,PlotRange->{{0,2},{0,0.05}},ImageSize->400]

k = 0.172

solution 1 (小)

(*---Part 1:Re-solve the system with the found parameters---*)
ClearAll["Global`*"];
(*Extract the final values for rir and a0 from the FindRoot result*)
finalrir=0.0202116;
finala0=0.000511577;
α:=0.1
k:=0.172
uv:=100
color={Yellow,Dashed};
lg={"k="<>ToString[k]}
(*B:=0.06*)
(*Solve the ODEs one last time with the optimal parameters to get the solution functions*)
finalSolution=NDSolve[{L''[ρ]==(-2 L[ρ] (2 k^2 α a[ρ]^2+2 α ρ^4 a'[ρ]^2)-3 ρ^5 L'[ρ])/ρ^6,a''[ρ]==-((-ρ^4 (k^2 a[ρ] (ρ^2-4 α L[ρ]^2)+ρ^3 (-3 ρ^2+4 α L[ρ]^2) a'[ρ])-8 α ρ^8 L[ρ] a'[ρ] L'[ρ])/(ρ^8 (ρ^2-4 α L[ρ]^2))),L[finalrir]==finalrir,L'[finalrir]==0,a[finalrir]==finala0,a'[finalrir]==0},{L,a},{ρ,finalrir,uv},Method->"StiffnessSwitching"];

(*Get the final interpolating functions*)
LsolFinal=L/. finalSolution[[1]];
asolFinal=a/. finalSolution[[1]];

(*Plot L(ρ) and k*a(ρ) on the same graph for comparison*)
plotL1=Plot[LsolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","L[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color,PlotLegends->lg]
plota1=Plot[k*asolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","k a[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color,PlotLegends->lg]
(*Plot the expression k^2*a(ρ)^2/ρ^3*)
plot2=Plot[-4α((((k^2)*((asolFinal[ρ])^2))/(ρ^4))+((asolFinal'[ρ])^2)),{ρ,finalrir,1},PlotStyle->color,PlotLegends->lg,AxesLabel->{"ρ","-4α((((k^2)*((a[ρ])^2))/(ρ^4))+a'[ρ]^2)"},PlotRange->{{0,1},Full}]

solution  2

(*(*---Part 1:Re-solve the system with the found parameters---*)
ClearAll["Global`*"];
*)
color=Green
(*Extract the final values for rir and a0 from the FindRoot result*)
finalrir=0.04802769164300066`;
finala0=0.041687368293961274`;
α:=0.1
k:=0.172
uv:=100
color={Yellow};
(*B:=0.06*)
(*Solve the ODEs one last time with the optimal parameters to get the solution functions*)
finalSolution=NDSolve[{L''[ρ]==(-2 L[ρ] (2 k^2 α a[ρ]^2+2 α ρ^4 a'[ρ]^2)-3 ρ^5 L'[ρ])/ρ^6,a''[ρ]==-((-ρ^4 (k^2 a[ρ] (ρ^2-4 α L[ρ]^2)+ρ^3 (-3 ρ^2+4 α L[ρ]^2) a'[ρ])-8 α ρ^8 L[ρ] a'[ρ] L'[ρ])/(ρ^8 (ρ^2-4 α L[ρ]^2))),L[finalrir]==finalrir,L'[finalrir]==0,a[finalrir]==finala0,a'[finalrir]==0},{L,a},{ρ,finalrir,uv},Method->"StiffnessSwitching"];

(*Get the final interpolating functions*)
LsolFinal=L/. finalSolution[[1]];
asolFinal=a/. finalSolution[[1]];

(*Plot L(ρ) and k*a(ρ) on the same graph for comparison*)
plotL2=Plot[LsolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","L[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color]
plota2=Plot[k*asolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","k a[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color]
(*Plot the expression k^2*a(ρ)^2/ρ^3*)
plot2=Plot[-4α((((k^2)*((asolFinal[ρ])^2))/(ρ^4))+((asolFinal'[ρ])^2)),{ρ,finalrir,1},PlotStyle->color,PlotLegends->{"k=0.172"},AxesLabel->{"ρ","-4α((((k^2)*((a[ρ])^2))/(ρ^4))+a'[ρ]^2)"},PlotRange->{{0,1},Full}]

Show[Plot[x,{x,0,0.06},PlotStyle->Black],plotL1,plotL2,PlotRange->{{0,1},{0,0.043}},ImageSize->400]
Show[plota1,plota2,PlotRange->{{0,2},{0,0.05}},ImageSize->400]

k = 0.1715

solution 1 （小）

(*---Part 1:Re-solve the system with the found parameters---*)
ClearAll["Global`*"];

color=Green
(*Extract the final values for rir and a0 from the FindRoot result*)
finalrir=0.0129865152993;
finala0=5.987177496199467*10^(-6);
α:=0.1
k:=0.1715
uv:=100
color={Green,Dashed};
lg={"k="<>ToString[k]}
(*B:=0.06*)
(*Solve the ODEs one last time with the optimal parameters to get the solution functions*)
finalSolution=NDSolve[{L''[ρ]==(-2 L[ρ] (2 k^2 α a[ρ]^2+2 α ρ^4 a'[ρ]^2)-3 ρ^5 L'[ρ])/ρ^6,a''[ρ]==-((-ρ^4 (k^2 a[ρ] (ρ^2-4 α L[ρ]^2)+ρ^3 (-3 ρ^2+4 α L[ρ]^2) a'[ρ])-8 α ρ^8 L[ρ] a'[ρ] L'[ρ])/(ρ^8 (ρ^2-4 α L[ρ]^2))),L[finalrir]==finalrir,L'[finalrir]==0,a[finalrir]==finala0,a'[finalrir]==0},{L,a},{ρ,finalrir,uv},Method->"StiffnessSwitching",AccuracyGoal->20,PrecisionGoal->20];

(*Get the final interpolating functions*)
LsolFinal=L/. finalSolution[[1]];
asolFinal=a/. finalSolution[[1]];

(*Plot L(ρ) and k*a(ρ) on the same graph for comparison*)
plotL1=Plot[LsolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","L[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color,PlotLegends->lg]
plota1=Plot[k*asolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","k a[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color,PlotLegends->lg]
(*Plot the expression k^2*a(ρ)^2/ρ^3*)
plot2=Plot[-4α((((k^2)*((asolFinal[ρ])^2))/(ρ^4))+((asolFinal'[ρ])^2)),{ρ,finalrir,1},PlotStyle->color,PlotLegends->{"k=0.1715"},AxesLabel->{"ρ","-4α((((k^2)*((a[ρ])^2))/(ρ^4))+a'[ρ]^2)"},PlotRange->{{0,1},Full}]

solution  2 （大）

(*(*---Part 1:Re-solve the system with the found parameters---*)
ClearAll["Global`*"];*)
(*Extract the final values for rir and a0 from the FindRoot result*)
finalrir=0.049414542328862436`;
finala0=0.045741162400818484`;
α:=0.1
k:=0.1715
uv:=100
color={Green};
(*B:=0.06*)
(*Solve the ODEs one last time with the optimal parameters to get the solution functions*)
finalSolution=NDSolve[{L''[ρ]==(-2 L[ρ] (2 k^2 α a[ρ]^2+2 α ρ^4 a'[ρ]^2)-3 ρ^5 L'[ρ])/ρ^6,a''[ρ]==-((-ρ^4 (k^2 a[ρ] (ρ^2-4 α L[ρ]^2)+ρ^3 (-3 ρ^2+4 α L[ρ]^2) a'[ρ])-8 α ρ^8 L[ρ] a'[ρ] L'[ρ])/(ρ^8 (ρ^2-4 α L[ρ]^2))),L[finalrir]==finalrir,L'[finalrir]==0,a[finalrir]==finala0,a'[finalrir]==0},{L,a},{ρ,finalrir,uv},Method->"StiffnessSwitching"];

(*Get the final interpolating functions*)
LsolFinal=L/. finalSolution[[1]];
asolFinal=a/. finalSolution[[1]];

(*Plot L(ρ) and k*a(ρ) on the same graph for comparison*)
plotL2=Plot[LsolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","L[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color]
plota2=Plot[k*asolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","k a[ρ]"},PlotRange->{{0,10},Full},PlotStyle->color]
(*Plot the expression k^2*a(ρ)^2/ρ^3*)
plot2=Plot[-4α((((k^2)*((asolFinal[ρ])^2))/(ρ^4))+((asolFinal'[ρ])^2)),{ρ,finalrir,uv},PlotStyle->color,PlotLegends->{"k=0.1715"},AxesLabel->{"ρ","-4α((((k^2)*((a[ρ])^2))/(ρ^4))+a'[ρ]^2)"},PlotRange->{{0,1},Full}]

Show[Plot[x,{x,0,0.06},PlotStyle->Black],plotL1,plotL2,PlotRange->{{0,1},{0,0.043}},ImageSize->400]
Show[plota1,plota2,PlotRange->{{0,2},{0,0.05}},ImageSize->400]

k=0 (Bz=0.05)

at k= 0, the eoms reduced to

Clear[k,α,Bz]
{L'[ρ] (-4 (Bz^2) β-3 (ρ^4)-4 (k^2) β ((a[ρ])^2)+8 (k^2) β ρ a[ρ] a'[ρ]+4 β (ρ^4) a'[ρ] (3 a'[ρ]+2 ρ a''[ρ]))+ρ (4 (Bz^2) β-(ρ^4)+4 (k^2) β ((a[ρ])^2)+4 β (ρ^4) ((a'[ρ])^2)) L''[ρ]==((2 L[ρ] (2 (Bz^2) α-γ (ρ^2)+2 (k^2) α ((a[ρ])^2)+2 α (ρ^4) ((a'[ρ])^2)))/ρ),(((4 α ((L[ρ])^2)+(ρ^2) (-1+4 β ((L'[ρ])^2))) ((k^2) a[ρ]-(ρ^4) a''[ρ]))/(ρ^3))==a'[ρ] (4 α ((L[ρ])^2)+8 α ρ L[ρ] L'[ρ]+(ρ^2) (-3+4 β L'[ρ] (3 L'[ρ]+2 ρ L''[ρ])))}/.{(*Bz->0,*)β->0,γ->0,k->0}
Solve[%,{L''[ρ],a''[ρ]}][[1]]

L''[ρ]->((-4 (Bz^2) α L[ρ]-4 α (ρ^4) L[ρ] ((a'[ρ])^2)-3 (ρ^5) L'[ρ])/(ρ^6))/.{a'[ρ]->0}

(*---Part 1:Re-solve the system with the found parameters---*)
ClearAll["Global`*"];
finalrir=0.10032867197814097`;
α:=0.1
Bz:=0.05
uv:=100
finalSolution=NDSolve[{-3 (ρ^4) L'[ρ]-(ρ^5) L''[ρ]==((2 L[ρ] (2 (Bz^2) α))/ρ),L[finalrir]==finalrir,L'[finalrir]==0},{L,a},{ρ,finalrir,uv}];

(*Get the final interpolating functions*)
LsolFinal=L/. finalSolution[[1]];
(*asolFinal=a/. finalSolution[[1]];*)

(*Plot L(ρ) and k*a(ρ) on the same graph for comparison*)
plot1=Plot[LsolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","L[ρ]"},PlotRange->{{0,10},Full},PlotStyle->Gray,PlotLegends->{"k=0"}]
(*plota=Plot[asolFinal[ρ],{ρ,finalrir,uv},AxesLabel->{"ρ","k a[ρ]"}]*)

all in one plot (figure 2a/2b)

Show all the output together:

Show[%408,%487,%215,%395,%304,%274,%131,Plot[x,{x,0,0.11},PlotStyle->Black],PlotRange->{{0,2},{0,0.11}}]
Show[%488,%216,%396,%305,%275,%132,PlotRange->{{0,2},{0,0.05}}]

Figure 3

EoMs

ClearAll["Global`*"];

(*equations of motion for L(ρ) and a(ρ) derived by the action eq.(4.1) in the letter*)eomL1=L''[ρ]==-(((-2 L[ρ] (2 (Bz^2) α-γ (ρ^2)+2 (k^2) α ((a[ρ])^2)+2 α (ρ^4) ((a'[ρ])^2))+ρ (-4 (Bz^2) β-3 (ρ^4)-4 (k^2) β ((a[ρ])^2)+8 (k^2) β ρ a[ρ] a'[ρ]+12 β (ρ^4) ((a'[ρ])^2)) L'[ρ]) (4 α ((L[ρ])^2)+(ρ^2) (-1+4 β ((L'[ρ])^2)))+8 β (ρ^2) a'[ρ] L'[ρ] ((k^2) a[ρ] (4 α ((L[ρ])^2)+(ρ^2) (-1+4 β ((L'[ρ])^2)))-(ρ^3) a'[ρ] (4 α ((L[ρ])^2)+8 α ρ L[ρ] L'[ρ]+3 (ρ^2) (-1+4 β ((L'[ρ])^2)))))/(-64 (β^2) (ρ^8) ((a'[ρ])^2) ((L'[ρ])^2)+(ρ^2) (4 (Bz^2) β-(ρ^4)+4 (k^2) β ((a[ρ])^2)+4 β (ρ^4) ((a'[ρ])^2)) (4 α ((L[ρ])^2)+(ρ^2) (-1+4 β ((L'[ρ])^2)))))
eoma1=a''[ρ]==(((k^2) a[ρ] ((ρ^2)-4 α ((L[ρ])^2))+(ρ^3) (-3 (ρ^2)+4 α ((L[ρ])^2)) a'[ρ]) (4 (Bz^2) β-(ρ^4)+4 (k^2) β ((a[ρ])^2)+4 β (ρ^4) ((a'[ρ])^2))+8 (ρ^4) L[ρ] a'[ρ] (8 (Bz^2) α β-2 β γ (ρ^2)-α (ρ^4)+8 α β ((k^2) ((a[ρ])^2)+(ρ^4) ((a'[ρ])^2))) L'[ρ]+4 β (ρ^2) ((k^2) a[ρ] (-4 (Bz^2) β+(ρ^4)-4 (k^2) β ((a[ρ])^2))+(ρ^3) (20 (Bz^2) β+3 (ρ^4)+20 (k^2) β ((a[ρ])^2)) a'[ρ]-20 (k^2) β (ρ^4) a[ρ] ((a'[ρ])^2)-12 β (ρ^7) ((a'[ρ])^3)) ((L'[ρ])^2))/((ρ^4) ((ρ^2)-4 α ((L[ρ])^2)) (4 (Bz^2) β-(ρ^4)+4 (k^2) β ((a[ρ])^2)+4 β (ρ^4) ((a'[ρ])^2))+4 β (ρ^6) (-4 (Bz^2) β+(ρ^4)-4 (k^2) β ((a[ρ])^2)+12 β (ρ^4) ((a'[ρ])^2)) ((L'[ρ])^2))

a''[ρ]==(((k^2) a[ρ] ((ρ^2)-4 α ((L[ρ])^2))+(ρ^3) (-3 (ρ^2)+4 α ((L[ρ])^2)) a'[ρ]) (4 (Bz^2) β-(ρ^4)+4 (k^2) β ((a[ρ])^2)+4 β (ρ^4) ((a'[ρ])^2))+8 (ρ^4) L[ρ] a'[ρ] (8 (Bz^2) α β-2 β γ (ρ^2)-α (ρ^4)+8 α β ((k^2) ((a[ρ])^2)+(ρ^4) ((a'[ρ])^2))) L'[ρ]+4 β (ρ^2) ((k^2) a[ρ] (-4 (Bz^2) β+(ρ^4)-4 (k^2) β ((a[ρ])^2))+(ρ^3) (20 (Bz^2) β+3 (ρ^4)+20 (k^2) β ((a[ρ])^2)) a'[ρ]-20 (k^2) β (ρ^4) a[ρ] ((a'[ρ])^2)-12 β (ρ^7) ((a'[ρ])^3)) ((L'[ρ])^2))/((ρ^4) ((ρ^2)-4 α ((L[ρ])^2)) (4 (Bz^2) β-(ρ^4)+4 (k^2) β ((a[ρ])^2)+4 β (ρ^4) ((a'[ρ])^2))+4 β (ρ^6) (-4 (Bz^2) β+(ρ^4)-4 (k^2) β ((a[ρ])^2)+12 β (ρ^4) ((a'[ρ])^2)) ((L'[ρ])^2))

Bz=1, B=1, γ = 0, α =  0.1, β = 0,  compare 1. a(r)=0 ; 2. B=1, k->0 or Infinity + Condensate, uv=100

Clear[α,β,γ,k,Bz,Et,B]

sub={Bz->1,α->0.1,γ->0,β->0};
eomL=eomL1/.sub//FullSimplify
eoma=eoma1/.sub//FullSimplify

(*table for the IR cutoff and a0 for varied k*)
({{"k", "ir",        "a0"}})
({{0.01`, 0.5335608579123652`, 99.99212398538721`}})
({{0.2`, 0.5299707353060042`, 4.843998017245209`}})
({{0.4`, 0.5197505409785911`, 2.1976916097333357`}})
({{0.7`, 0.49660693051533467`, 0.9516899044379133`}})
({{1, 0.4744376067105367`, 0.4359529593771699`}})
({{1.2`, 0.4643452611661266`, 0.25931381895754746`}})
({{1.6`, 0.4541515551692595`, 0.09344829721279928`}})
({{2, 0.4506528924564342`, 0.03470269449321098`}})
({{2.5`, 0.4492966817455199`, 0.010355445462531141`}})
({{5, 0.44869000751342897`, 0.0000294261107751334`}})
({{6, 0.4486826815921271`, 2.9273015665581446`*^-6}})
({{14, 0.4486789690723186`, 3.6385436151263124`*^-14}})

(*Bz:=1*)
k:=0.2
B:=1
uv:=100
guessIR:=0.5299707353060042`
guessa0:=1/k
(*Plot the final solution*)
color=Table[ColorData["Rainbow"][i],{i,Subdivide[6-1]}]
plotstyle=color[[1]];
legends={"k="<>ToString[k](*<>",B="<>ToString[B]<>",(Subscript[B, z])="<>ToString[Bz]*)}
(*Function to solve the ODEs and return the residuals*)
SolveWithParameters[ir_?NumericQ,a0_?NumericQ]:=Module[{solution,Lsol,asol},(*Define the initial conditions with the input parameters*)icsL={L[ir]==ir,L'[ir]==0};
icsa={a[ir]==a0,a'[ir]==0};
(*Solve the system of ODEs*)solution=NDSolve[{eomL,eoma,icsL,icsa},{L,a},{ρ,ir,uv},Method->"StiffnessSwitching",PrecisionGoal->20,AccuracyGoal->20];
(*Extract the functions from the solution*)Lsol=L/. solution[[1]];
asol=a/. solution[[1]];
(*Return the residuals we want to be zero*){asol[uv]*k-B,Lsol[uv]}];

(*Find the root of the function*)(*The first list {ir,0.5} gives the variable and an initial guess.*)(*The second list {a0,1.0} does the same for the other variable.*)foundParameters=FindRoot[SolveWithParameters[rir,a0],{{rir,guessIR},{a0,guessa0}}]
TableForm[{{k,foundParameters[[1,2]],foundParameters[[2,2]]}}]

(*Extract the final values for ir and a0 from the solution*)
finalIr=rir/. foundParameters;
finalA0=a0/. foundParameters;
(*---Final Solution with the correct initial conditions---*)
finalIcsL={L[finalIr]==finalIr,L'[finalIr]==0};
finalIcsa={a[finalIr]==finalA0,a'[finalIr]==0};

finalSolution=NDSolve[{eomL,eoma,finalIcsL,finalIcsa},{L,a},{ρ,finalIr,uv}];

(*Check the final values to confirm they are correct*)
Print["Final value of a[uv]*k: ",(a[uv]*k)/. finalSolution];
Print["Final value of L[uv]: ",L[uv]/. finalSolution];

Plot[Evaluate[L[ρ]/. finalSolution],{ρ,finalIr,uv},PlotLegends->legends,PlotRange->Full,AxesLabel->{"ρ","L(ρ)"},PlotStyle->plotstyle]

(*Plot the final solution*)
Plot[Evaluate[k*a[ρ]/. finalSolution],{ρ,finalIr,uv},PlotLegends->legends,PlotRange->All,AxesLabel->{"ρ","k*a(ρ)"},PlotStyle->plotstyle]

(*Get the final interpolating functions*)
α:=0.1
Bz:=1
LsolFinal=L/. finalSolution[[1]];
asolFinal=a/. finalSolution[[1]];
plot2=Plot[-4α((((k^2)*((asolFinal[ρ])^2))/(ρ^4))+((asolFinal'[ρ])^2)),{ρ,finalIr,15},PlotStyle->plotstyle,PlotLegends->{"k="<>ToString[k]},AxesLabel->{"ρ","-4α((((k^2)*((a[ρ])^2))/(ρ^4))+a'[ρ]^2)"},PlotRange->{{0,15},Full}]
plot2=Plot[-4α(((Bz^2)/(ρ^4))+(((k^2)*((asolFinal[ρ])^2))/(ρ^4))+((asolFinal'[ρ])^2)),{ρ,finalIr,15},PlotStyle->plotstyle,PlotLegends->{"k="<>ToString[k]},AxesLabel->{"ρ","-4α(((Bz^2)/(ρ^4))+(((k^2)*((a[ρ])^2))/(ρ^4))+a'[ρ]^2)"},PlotRange->{{0,15},Full}]

(*plots is calclated by different k, with given IR cut off and a0 in the table above*)
Show[Plot[x,{x,0,0.6},PlotStyle->Black],%803,%909,%1517,%1781,%1424,PlotRange->{{0,3},{0,0.6}}]
Show[%804,%910,%1518,%1782,%1425,PlotRange->{{0,10},{0,1}},AxesOrigin->{0,0}]
Show[%814,%915,%1523,%1793,%1442]
Show[%815,%916,%1524,%1794,%1443,PlotRange->{{0,4},Full}]

Show[%1922,AxesLabel->{"ρ",None}](*(mL^2)-(-((ΔmL)^2))*)

Show[%1920,AxesLabel->{"ρ","L"}]

(Legended[, ({{(LineLegend["k=0.2"])}, {(LineLegend["k=0.7"])}, {(LineLegend["k=1.6"])}, {(LineLegend["k=6"])}, {(LineLegend["k=14"])}})])

Show[%147,%100,PlotRange->{{0,3},Full}]

condensate

Bz:=1;
α:=0.1;
uv:=100;
B:=1;
(*value of k*)
kList={0.01,0.1,0.5,1,2,4,6,8(*,14*)};

(*save the results*)
results=<||>;
plotDataL={};
plotDataA={};
(*---New---Initialize a list to store the results for the final plot*)
calculatedData={};


getInitialGuess[k_]:=Module[{prevK,prevIr,prevA0,prevProduct},If[Length[Keys[results]]>0,prevK=Last[Keys[results]];
{prevIr,prevA0}=Lookup[results[prevK],{"ir","a0"}];
prevProduct=prevK*prevA0;
Which[k<0.1,{0.53,1/k},k<1,{0.53,prevProduct/k},True,{0.45,prevProduct/k}],If[k<0.1,{0.53,1/k},{0.45,0.1/k}]]];


Do[Print["========================================"];
Print["Solving for k = ",k];
{currentIrGuess,currentA0Guess}=getInitialGuess[k];
Print["Initial guess: ir = ",currentIrGuess,", a0 = ",currentA0Guess];
SolveWithParameters[ir_?NumericQ,a0_?NumericQ]:=Module[{solution,Lsol,asol},icsL={L[ir]==ir,L'[ir]==0};
icsa={a[ir]==a0,a'[ir]==0};
solution=NDSolve[{eomL,eoma,icsL,icsa},{L,a},{ρ,ir,uv},Method->"StiffnessSwitching"];
Lsol=L/. solution[[1]];
asol=a/. solution[[1]];
{asol[uv]*k-B,Lsol[uv]}];
foundParameters=Quiet@Check[FindRoot[SolveWithParameters[rir,a0],{{rir,currentIrGuess,0.5353644798214199`,19},{a0,currentA0Guess,0,Max[10,10/k]}},MaxIterations->50],$Failed];
If[foundParameters===$Failed,Print["FindRoot failed for k = ",k,". Trying fallback method..."];
foundParameters=Quiet@Check[FindRoot[SolveWithParameters[rir,a0],{{rir,0.5,0.4,0.6},{a0,0.5/k,0,10/k}},MaxIterations->100],$Failed];];
If[foundParameters===$Failed,Print["Failed to find solution for k = ",k];
Continue[]];
finalIr=rir/. foundParameters;
finalA0=a0/. foundParameters;
finalIcsL={L[finalIr]==finalIr,L'[finalIr]==0};
finalIcsa={a[finalIr]==finalA0,a'[finalIr]==0};
finalSolution=NDSolve[{eomL,eoma,finalIcsL,finalIcsa},{L,a},{ρ,finalIr,uv}];
(*---New Calculation Section---*)(*Check if the solution was found before proceeding*)If[finalSolution=!={},(*Extract the solution for L from the numerical solution*)Lsol=L/. finalSolution[[1]];
(*Evaluate the derivative L' at the boundary rho=uv*)LprimeAtUv=Lsol'[uv];
(*Calculate the desired quantity*)calculatedValue=-uv^3*LprimeAtUv;
(*Add the {k,calculatedValue} pair to our data list*)AppendTo[calculatedData,{k,calculatedValue}];];
(*---End of New Calculation Section---*)results[k]=<|"ir"->finalIr,"a0"->finalA0,"k*a0"->k*finalA0,"L"->(L/. finalSolution[[1]]),"a"->(a/. finalSolution[[1]])|>;
AppendTo[plotDataL,{k,L/. finalSolution[[1]]}];
AppendTo[plotDataA,{k,a/. finalSolution[[1]]}];
Print["Solved: k = ",k,", ir = ",finalIr,", a0 = ",finalA0,", k*a0 = ",k*finalA0],{k,kList}];


(*plot L(ρ)*)
allLPlots=Show[Table[Plot[Evaluate[results[k]["L"][ρ]],{ρ,results[k]["ir"],uv},PlotStyle->ColorData[97][Position[kList,k][[1,1]]],PlotRange->{{0,5},Full}],{k,kList}],AxesLabel->{"ρ","L(ρ)"},ImageSize->600,PlotLegends->LineLegend[kList,LegendLabel->"k values"],PlotRange->{{0,10},Full}];

(*plot k*a(ρ)*)
allAPlots=Show[Table[Plot[Evaluate[k*results[k]["a"][ρ]],{ρ,results[k]["ir"],uv},PlotStyle->ColorData[97][Position[kList,k][[1,1]]],PlotRange->{{0,uv},Full}],{k,kList}],AxesLabel->{"ρ","k*a(ρ)"},PlotLabel->"k*a(ρ) for various k values",ImageSize->600,PlotLegends->LineLegend[kList,LegendLabel->"k values"],PlotRange->{{0,uv},Full}(*Let Show determine the overall range based on individual plots*)];

(*plot*)
parameterPlot=Grid[{{ListPlot[{#1,#2}&@@@(Thread[{kList,results[#]["ir"]&/@kList}]),Joined->True,PlotRange->All,ImageSize->300,AxesLabel->{"k","ir"},PlotLabel->"ir vs k"],ListPlot[{#1,#2}&@@@(Thread[{kList,results[#]["a0"]&/@kList}]),Joined->True,PlotRange->All,ImageSize->300,ScalingFunctions->"Log",AxesLabel->{"k","a0"},PlotLabel->"a0 vs k (log scale)"]},{ListPlot[{#1,#2}&@@@(Thread[{kList,results[#]["k*a0"]&/@kList}]),Joined->True,PlotRange->All,ImageSize->300,AxesLabel->{"k","k*a0"},PlotLabel->"k*a0 vs k"],ListPlot[{#1,#2}&@@@(Thread[{kList,results[#]["a0"]&/@kList}]),Joined->True,PlotRange->All,ImageSize->300,ScalingFunctions->{"Log","Log"},AxesLabel->{"k","a0"},PlotLabel->"a0 vs k (log-log scale)"]}}];

parameterTable=Table[{k,results[k]["ir"],results[k]["a0"],results[k]["k*a0"]},{k,kList}];
table=TableForm[parameterTable,TableHeadings->{None,{"k","ir","a0","k*a0"}}];

(*---New---Plot the final calculated data*)
finalPlot=ListLinePlot[calculatedData,PlotRange->All,AxesLabel->{"k","⟨(OverBar[q])q⟩"},PlotMarkers->{Automatic,Medium},GridLines->Automatic(*,PlotLabel->"k vs. Calculated Boundary Condition"*),ImageSize->Large]

(*Print[table];
Print[allLPlots];
Print[allAPlots];*)
Print[finalPlot];

k->0 limit+condensate

sub={(*α->0.1,*)γ->0,β->0};
eomL=eomL1/.sub//FullSimplify
eoma=eoma1/.sub//FullSimplify

ClearAll["Global`*"]
L''[ρ]==((-4 α L[ρ] ((Bz^2)+(k^2) ((a[ρ])^2)+(ρ^4) ((a'[ρ])^2))-3 (ρ^5) L'[ρ])/(ρ^6))/.{a[ρ]->B/k,a'[ρ]->0}
%/.B->1

(*Define the parameters*)
ir:=0.5335699688080151`;
uv:=100;
Bz:=1;
α:=0.1;
color={Black,Dashed};
icsL={L[ir]==ir,L'[ir]==0};
solution1=NDSolve[{L''[ρ]==((-4 (1+(Bz^2)) α L[ρ]-3 (ρ^5) L'[ρ])/(ρ^6)),icsL},L,{ρ,ir,uv},Method->"StiffnessSwitching"];
Plot[Evaluate[L[ρ]/. solution1],{ρ,ir,uv},AxesLabel->{"ρ","L(ρ)"},PlotLegends->{"k=0"(*"(Subscript[B, z])=1,a(ρ)=0"*)},PlotRange->All,PlotStyle->color]

Figure 4

EoMs

ClearAll["Global`*"];

(*eomL1 and eoma1 are Euler Lagrangian equation of the action i.e. equation (4.1) in the letter*)
eomL1=L''[ρ]==-(((-2 L[ρ] (2 (Bz^2) α-γ (ρ^2)+2 (k^2) α ((a[ρ])^2)+2 α (ρ^4) ((a'[ρ])^2))+ρ (-4 (Bz^2) β-3 (ρ^4)-4 (k^2) β ((a[ρ])^2)+8 (k^2) β ρ a[ρ] a'[ρ]+12 β (ρ^4) ((a'[ρ])^2)) L'[ρ]) (4 α ((L[ρ])^2)+(ρ^2) (-1+4 β ((L'[ρ])^2)))+8 β (ρ^2) a'[ρ] L'[ρ] ((k^2) a[ρ] (4 α ((L[ρ])^2)+(ρ^2) (-1+4 β ((L'[ρ])^2)))-(ρ^3) a'[ρ] (4 α ((L[ρ])^2)+8 α ρ L[ρ] L'[ρ]+3 (ρ^2) (-1+4 β ((L'[ρ])^2)))))/(-64 (β^2) (ρ^8) ((a'[ρ])^2) ((L'[ρ])^2)+(ρ^2) (4 (Bz^2) β-(ρ^4)+4 (k^2) β ((a[ρ])^2)+4 β (ρ^4) ((a'[ρ])^2)) (4 α ((L[ρ])^2)+(ρ^2) (-1+4 β ((L'[ρ])^2)))))(*copy from eomL1*)
eoma1=a''[ρ]==(((k^2) a[ρ] ((ρ^2)-4 α ((L[ρ])^2))+(ρ^3) (-3 (ρ^2)+4 α ((L[ρ])^2)) a'[ρ]) (4 (Bz^2) β-(ρ^4)+4 (k^2) β ((a[ρ])^2)+4 β (ρ^4) ((a'[ρ])^2))+8 (ρ^4) L[ρ] a'[ρ] (8 (Bz^2) α β-2 β γ (ρ^2)-α (ρ^4)+8 α β ((k^2) ((a[ρ])^2)+(ρ^4) ((a'[ρ])^2))) L'[ρ]+4 β (ρ^2) ((k^2) a[ρ] (-4 (Bz^2) β+(ρ^4)-4 (k^2) β ((a[ρ])^2))+(ρ^3) (20 (Bz^2) β+3 (ρ^4)+20 (k^2) β ((a[ρ])^2)) a'[ρ]-20 (k^2) β (ρ^4) a[ρ] ((a'[ρ])^2)-12 β (ρ^7) ((a'[ρ])^3)) ((L'[ρ])^2))/((ρ^4) ((ρ^2)-4 α ((L[ρ])^2)) (4 (Bz^2) β-(ρ^4)+4 (k^2) β ((a[ρ])^2)+4 β (ρ^4) ((a'[ρ])^2))+4 β (ρ^6) (-4 (Bz^2) β+(ρ^4)-4 (k^2) β ((a[ρ])^2)+12 β (ρ^4) ((a'[ρ])^2)) ((L'[ρ])^2))(*copy from eoma1*)

γ=-1.7, k≠0, α=0, β= + 0.1, Bz=0,

Clear[α,β,γ,k,Bz,Et]
sub={γ->-1.7,α->0(*,k->0*),Bz->0,Et->0}
eomL=eomL1/.sub//FullSimplify
eoma=eoma1/.sub//FullSimplify

k=0 limit

(*Define constants*)
Clear[γ,B,β]
γ=-1.7;
B=1;
β=0.1;

(*Define a function that returns the value of the solution at rho=100*)
endValue[ir_?NumericQ]:=NDSolveValue[{((2 γ ρ L[ρ]-(3 (ρ^4)+4  β (B^2)) L'[ρ])/((ρ^5)-4  β ρ (B^2)))==L''[ρ],L[ir]==ir,L'[ir]==0},L[100],(*This is the value we want to be zero*){ρ,ir,100},Method->"StiffnessSwitching"];
(*Find the value of ir that makes endValue[ir]==0*)irSolution=FindRoot[endValue[x]==0,{x,1}]

(*γ:=-1.7;
B:=1;
β:=0.1;*)
ir=0.8084574296049702`
Lsol=NDSolveValue[{((2 γ ρ L[ρ]-(3 (ρ^4)+4  β (B^2)) L'[ρ])/((ρ^5)-4  β ρ (B^2)))==L''[ρ],L[ir]==ir,L'[ir]==0},L,{ρ,ir,100}];
Plot[Lsol[x],{x,ir,100},PlotStyle->{Dashed,Black},PlotRange->{{0,4},Full},PlotLegends->{"k=0"}]

eom for k->0

Clear[α,β,γ,k,Bz,Et,B]
sub={(*γ->-1.7,*)α->0(*,k->0*),Bz->0,Et->0}
{L'[ρ] (-4 (Bz^2) β-3 (ρ^4)-4 (k^2) β ((a[ρ])^2)+8 (k^2) β ρ a[ρ] a'[ρ]+4 β (ρ^4) a'[ρ] (3 a'[ρ]+2 ρ a''[ρ]))+ρ (4 (Bz^2) β-(ρ^4)+4 (k^2) β ((a[ρ])^2)+4 β (ρ^4) ((a'[ρ])^2)) L''[ρ]==((2 L[ρ] (2 (Bz^2) α-γ (ρ^2)+2 (k^2) α ((a[ρ])^2)+2 α (ρ^4) ((a'[ρ])^2)))/ρ),(((4 α ((L[ρ])^2)+(ρ^2) (-1+4 β ((L'[ρ])^2))) ((k^2) a[ρ]-(ρ^4) a''[ρ]))/(ρ^3))==a'[ρ] (4 α ((L[ρ])^2)+8 α ρ L[ρ] L'[ρ]+(ρ^2) (-3+4 β L'[ρ] (3 L'[ρ]+2 ρ L''[ρ])))}/.sub

finetune ir

(*this is the table of the finetuned value of IR cutoff with varied k*)
(*k,    the IR cutoff of ρ,   a0=a[IR]*)
({{0.1`, 0.8081733441207698`, 9.966010783898607`}})
({{0.5`, 0.8022659222869759`, 1.8320270382006123`}})
({{2, 0.7750944269371997`, 0.14573214840756846`}})
({{2.5`, 0.7715043066541313`, 0.06872729316202988`}})
({{3, 0.769459862928202`, 0.032864492788946195`}})
({{3.5`, 0.7683039070136229`, 0.015893257283430474`}})
({{5, 0.7670432135134589`, 0.0018859445149311735`}})
({{15, 0.7666712264637391`, 2.4119321572555244`*^-9}})

(*Define the parameters*)
Clear[α,β,k,ir,uv]
B:=1(*K*(Subscript[a, uv]) so that to fix the boundary magnetic field*)
uv:=100;(*uv cutoff of ρ*)
β:=0.1;
k:=5;
irg:=0.809;(*use the IR cutoff of ρ*)
a0g:=(*1/k*)4.274*E^(-1.785*k); (*for large k use 4.274*E^(-1.785*k) to guess a0; for small k use 1/k to guess the a0, or use the values from the above table*)

color=Table[ColorData["Rainbow"][i],{i,Subdivide[9-1]}];
plotstyle=color[[6]];(*ser the color of the plot by hand*)

(*Function to solve the ODE*)
SolveWithParameters[ir_?NumericQ,a0_?NumericQ]:=Module[{solution,Lsol,asol},(*Define the initial conditions with the input parameters*)icsL={L[ir]==ir,L'[ir]==0};
icsa={a[ir]==a0,a'[ir]==0};
(*Solve the system of ODEs*)solution=NDSolve[{eomL,eoma,icsL,icsa},{L,a},{ρ,ir,uv},Method->"StiffnessSwitching",PrecisionGoal->20,AccuracyGoal->20];
(*Extract the functions from the solution*)Lsol=L/. solution[[1]];
asol=a/. solution[[1]];
(*Return the residuals we want to be zero*){asol[uv]*k-B,Lsol[uv]}];
foundParameters=FindRoot[SolveWithParameters[rir,a0],{{rir,irg},{a0,a0g}}]
TableForm[{{k,foundParameters[[1,2]],foundParameters[[2,2]]}}]
(*Extract the final values for ir and a0 from the solution*)finalIr=rir/. foundParameters;
finalA0=a0/. foundParameters;

(*---Final Solution with the correct initial conditions---*)
finalIcsL={L[finalIr]==finalIr,L'[finalIr]==0};
finalIcsa={a[finalIr]==finalA0,a'[finalIr]==0};

finalSolution=NDSolve[{eomL,eoma,finalIcsL,finalIcsa},{L,a},{ρ,finalIr,uv}];

(*Check the final values to confirm they are correct*)
Print["Final value of a[uv]*k: ",(a[uv]*k)/. finalSolution];
Print["Final value of L[uv]: ",L[uv]/. finalSolution];

Plot[Evaluate[L[ρ]/. finalSolution],{ρ,finalIr,uv},PlotLegends->{"k="<>ToString[k]},PlotRange->Full,AxesLabel->{"ρ","L(ρ)"},PlotStyle->plotstyle]
(*Plot the final solution*)
Plot[Evaluate[k*a[ρ]/. finalSolution],{ρ,finalIr,uv},PlotLegends->{"k="<>ToString[k]},PlotRange->All,AxesLabel->{"ρ","a(ρ)"},PlotStyle->plotstyle]
Lsol=L/. finalSolution[[1]];
(*Evaluate the derivative L' at the boundary rho=uv*)LprimeAtUv=Lsol'[uv];
(*Calculate the desired quantity*)
Print[{"k","condensate"}](*return k and the condensate*)
{k,-uv^3*LprimeAtUv}

(*show all the plots with the finetuned value*)
Show[%225,%148,%92,%27,%16,Plot[x,{x,0,0.9},PlotStyle->Black],PlotRange->{{0,4},Full},AxesOrigin->{0,0}]
Show[%226,%149,%93,%28,PlotRange->{{0,20},Full},AxesOrigin->{0,0},AxesLabel->{"ρ","k*a(ρ)"}]

plot1=Show[,PlotRange->{{0.65,1},{0.65,0.82}},PlotLegends->False]
 plot2=Show[%343,%231,,PlotRange->{{0,4},{0,0.8}}]
insetPosition={3.7,0.8};
insetSize=1.5;
finalPlot=Show[plot2,Epilog->{Inset[plot1,insetPosition,{Right,Top},insetSize]},AxesLabel->{"ρ","L(ρ)"}]

qbq={{0.01`,1.808807347543084`},{0.5`,1.884744964703899`},{1,1.986135575765834`},{2,2.0717334125491607`},{2.5`,2.082617690038283`},{5,2.0869145436951295`},{10,2.086181106060439`}(*,{15,2.085999499577563`}*)};(*k against condensate*)
(*---New---Plot the final calculated data*)
ListLinePlot[qbq,PlotRange->All,PlotMarkers->{Automatic,Medium},GridLines->Automatic,ImageSize->400,AxesLabel->{"k","⟨(OverBar[q])q⟩"}]