ClearAll["Global'*"]

(*ResourceFunction["MaTeXInstall"][]*)
<<MaTeX`

A matching theory to characterize sound emission during vortex reconnection in quantum fluids: Supplemental material

The geometry of the reconnecting filaments

In this note we develop further the analytical calculations performed in Villois et al PRFluids 2017  [A. Villois, D. Proment, and G. Krstulovic, Phys. Rev. Fluids 2, 044701 (2017)] within the framework of the linear approximation.

The reconnecting wave-function of the filaments

Differently from Villois et al PRFluids 2017, we consider here a more complete ansatz for the wavefunction at the reconnection time (t=0). It reads

ψ0[x_,y_,z_]=(p ( z-(A/2 (Cos[θ] x+Sin[θ] y)^2+B/2 (-Sin[θ]x+Cos[θ]y)^2)/ζ)+I(z-(C/2 x^2+D/2 y^2)/ζ))/ζ^(5/2)
RestrictPar=p^2==1 && A∈Reals && B∈Reals && C∈Reals && D∈Reals && 0≤θ<π && ζ>0

Note that the parameters A, B, C, D θ, p, are non-dimensional numbers, and because it must be [|ψ|^2]=1/length^3, then [ζ]=length. This is the most general ansatz one can consider for reconnecting vortex filaments, at t=0, Taylor expanded up to the quadratic order in the spatial coordinates, where p=(-1, 1) controls the orientation of the real and imaginary parts.

Note also that the vortex filaments are the nodal lines of the field ψ0, that is the intersection of the paraboloid obtained by setting Re(ψ0)=0 and the hyperboloid obtained by setting Im(ψ0)=0. An example is shown below.

ParamValues={A->-2., B->-1., C->-2., D->1., θ->π/3., p->-1, ζ->1.};
RestrictPar/.ParamValues
Reψ0[x_, y_, z_]=((p (z-(((1/2) B (((y Cos[θ]-x Sin[θ]))^2)+(1/2) A (((x Cos[θ]+y Sin[θ]))^2))/ζ)))/(ζ^(5/2)))/.ParamValues;
Imψ0[x_, y_, z_]=(((z-((((C (x^2))/2)+((D (y^2))/2))/ζ)))/(ζ^(5/2)))/.ParamValues;
ContourPlot3D[{Reψ0[x, y, z]==0, Imψ0[x, y, z]==0}, {x,-2,2},{y,-2,2},{z,-2,2}, PlotLegends->{"Re(ψ0)==0", "Im(ψ0)==0"}, AxesLabel->Automatic, PlotLabel->ParamValues]
genIC=ContourPlot3D[{Reψ0[x, y, z]==0, Imψ0[x, y, z]==0}, {x,-2,2},{y,-2,2},{z,-2,2}, PlotLegends->{MaTeX[#, Magnification->1.5]&/@{"\\text{Re}(\\psi_\\text{r})=0","\\text{Im}(\\psi_\\text{r})=0"}}, AxesLabel->{MaTeX["x/\\zeta", Magnification->1.5], MaTeX["y/\\zeta", Magnification->1.5], MaTeX["z/\\zeta", Magnification->1.5]}, BaseStyle->{FontSize->15}];

The time evolution of the initial configuration about the reconnection is obtained by solving the (linear) Schrodinger equation. That implies ψ[t,x,y,z]=exp[i t (Γ/4π) Δ ]ψ0[x,y,z], where we assume Γ>0. As the wave function is quadratic at the reconnection time, it follows

ψ[t_,x_,y_,z_]=FullSimplify[Collect[Expand[ψ0[x,y,z]+ I t (Γ/(4π)) Laplacian[ψ0[x,y,z],{x,y,z}]], {I}]]
RestrictGamma=Γ>0

The vortex lines are found by finding the nodal lines of ψ. Taking real an imaginary part equal to 0, we obtain

Eq1=-B p π (x^2)-B p π (y^2)+D t Γ-A p (π ((x^2)+(y^2)))+C t Γ+4 p π z ζ+(A-B) p π ((-(x^2)+(y^2)) Cos[2 θ]-2 x y Sin[2 θ])//FullSimplify
Eq2=-2 D π (y^2)- B p t Γ-A p t Γ-C 2 π (x^2)+4 π z ζ//FullSimplify

These can be expressed in terms of z as

Collect[Solve[Eq1==0,z][[1]]//FullSimplify,{t}]
Collect[Solve[Eq2==0,z][[1]]//FullSimplify,{t}]

Note that the vortex filaments do not lay generally onto the x-y plane. The position in the z coordinate can be calculated for instance by the first equation

ReplaceZ=Collect[Solve[Eq1==0,z][[1]]//FullSimplify,{t}]

When projected onto the z=0 plane, the reconnecting vortex filaments follow the branches of a hyperbola, as already discussed already in Villois et al PRFluids 2017. Its equation reads

Hyperbola=Eq2/.ReplaceZ//Expand;
HyperbolaSimp=Collect[Hyperbola,{x^2,y^2,t Γ}]==0//FullSimplify

ParamValues={A->-2., B->-1, C->-2., D->1., θ->π/3., p->-1, ζ->1.};
RestrictPar/.ParamValues
hyp[x_, y_, t_]=Hyperbola/.ParamValues/.{Γ->1};
ContourPlot[{hyp[x, y, -1]==0, hyp[x, y, 1]==0}, {x, -2, 2}, {y, -2, 2}, PlotLegends->{"t=-1", "t=1"}, Axes->True, FrameLabel->Automatic, PlotLabel->ParamValues]
genICProj=ContourPlot[{hyp[x, y, -1]==0, hyp[x, y, 1]==0}, {x, -2, 2}, {y, -2, 2}, PlotLegends->{MaTeX["t=-1", Magnification->1.5], MaTeX["t=1", Magnification->1.5]}, Axes->True, FrameLabel->{MaTeX["x/\\zeta", Magnification->1.5], MaTeX["y/\\zeta", Magnification->1.5]}];

In particular, when θ=0, the equation of the hyperbola becomes more evident

Replaceθ0=θ->0;
HyperbolaSimp=Collect[HyperbolaSimp/.Replaceθ0, {x^2,y^2,t Γ}]//FullSimplify

Region of validity of the other free parameters

We denote R1m and R2m the vortices before the reconnection and R1p and R2p the vortices after the reconnection. Without loss of generality we assume that before the reconnection R1m ⊂ { y>0 } and R2m ⊂ { y<0 }, whereas after reconnection  R1p ⊂ { x>0 } and R2p ⊂ { x<0 }. Also, we assume that the vorticity orientation is such that the R1m points towards the negative x-direction while R2m towards the positive one; the orientations after follows reconnection rules so that R1p points towards the positive y-direction while R2p towards the negative one. For a visual representation see see Fig. 4 of the manuscript.

Condition for the existence of the hyperbola

At reconnection time, the projection of the vortex filaments onto the z=0 plane gives the hyperbola asymptotes

asympt=Solve[HyperbolaSimp/.{t->0}, {y}]

In order to have real asymptotes, the following restriction must be satisfied

RestrictAsympt=Reduce[{(((Sqrt[-A+C])/(Sqrt[B-D])))^2≥0}]//FullSimplify

Reduce[{RestrictAsympt, RestrictPar, RestrictGamma}]

Convention on the positions of the filaments

To fulfil the position convention R1m ⊂ { y>0 } and R2m ⊂ { y<0 }, the parameters must satisfy the restriction

Solve[HyperbolaSimp/.{x->0, t->-1}, y]

(((I (Sqrt[C+D+A (p^2)+B (p^2)]) (Sqrt[Γ]))/((Sqrt[p]) (Sqrt[2 B π-2 D π]))))^2>0//FullSimplify
RestrictFilBef=Reduce[{(((I (Sqrt[C+D+A (p^2)+B (p^2)]) (Sqrt[Γ]))/((Sqrt[p]) (Sqrt[2 B π-2 D π]))))^2>0}]//FullSimplify

To fulfil the position convention R1p ⊂ { x>0 } and R2p ⊂ { x<0 }, the parameters must satisfy the restriction

Solve[HyperbolaSimp/.{y->0, t->1}, x]

((((Sqrt[C+D+A (p^2)+B (p^2)]) (Sqrt[Γ]))/((Sqrt[p]) (Sqrt[2 A π-2 C π]))))^2>0//FullSimplify
RestrictFilAft=Reduce[{((((Sqrt[C+D+A (p^2)+B (p^2)]) (Sqrt[Γ]))/((Sqrt[p]) (Sqrt[2 A π-2 C π]))))^2>0}]//FullSimplify

Reduce[{RestrictFilAft, RestrictFilBef, RestrictAsympt, RestrictPar, RestrictGamma}]

Convention on the vorticity orientation of the filaments

One can use the pseudo-vorticity to estimate the vorticity direction of the filaments. The pseudo-vorticity is proportional to the cross production of the gradients of the real and imaginary parts of the wave function, resulting in

pseudoVort=Cross[Grad[Eq1, {x, y, z}], Grad[Eq2, {x, y, z}]]/.Replaceθ0//Expand
Collect[pseudoVort[[1]], {x, y}]
Collect[pseudoVort[[2]], {x, y}]

In order to have the filaments’ vorticity direction as the one sketched in Fig. 4 of the manuscript, the following condition must hold

{-16 B p (π^2) ζ+16 D p (π^2) ζ<0, 16 A p (π^2) ζ-16 C p (π^2) ζ<0}//FullSimplify
RestrictOr=Reduce[{-16 B p (π^2) ζ+16 D p (π^2) ζ<0, 16 A p (π^2) ζ-16 C p (π^2) ζ<0}//FullSimplify]

Reduce[{RestrictOr, RestrictFilAft, RestrictFilBef, RestrictAsympt, RestrictPar, RestrictGamma}]//FullSimplify

Combining all the restrictions

When combining all the restrictions one obtains

Restrictions[p_]=Reduce[{RestrictOr, RestrictFilAft, RestrictFilBef, RestrictAsympt, RestrictPar, RestrictGamma}]//FullSimplify

FullSimplify[Restrictions[-1]]
FullSimplify[Restrictions[1]]

The geometrical parameters of the filaments

We related here the set of the free parameters of the initial condition to the geometrical parameters of the filaments.

The reconnecting angle

Projected onto the z=0 plane,  the filaments form a hyperbola having asymptotes

asympt

This means that the reconnection angle after the reconnection, ϕ+, results in

ϕp=2 ArcTan[Sqrt[((-A+C)/(B-D))]]

Fixing the geometrical value of ϕ+Star results in the following constraint between the parameters

Solve[ϕp==ϕpStar, C]

ReplaceTan2Star={((Tan[((ϕpStar)/2)])^2)->Tan2Star};
ReplaceCStar={C->A+B ((Tan[((ϕpStar)/2)])^2)-D ((Tan[((ϕpStar)/2)])^2)}/.ReplaceTan2Star

The concavity

Let us now find out how the filaments are projected onto the  y=0 plane.

Solve[Eq1==0/.Replaceθ0/.{y^2->ySquare}, ySquare]
ReplaceYSquare={y^2->((-2 A p π (x^2)+C t Γ+D t Γ+4 p π z ζ)/(2 B p π))}

Solve[Eq2==0/.Replaceθ0/.ReplaceYSquare, z]

Let us now analyse the z-dependence of the filaments

Zvort=FullSimplify[((2 B C p π (x^2)-2 A D p π (x^2)+C D t Γ+(D^2) t Γ+A B (p^2) t Γ+(B^2) (p^2) t Γ)/(4 (B-D) p π ζ))]

Note that this is quadratic with respect to the x-directions, with concavity amplitude and sign controlled by the parameters A, B, C and D as

conc=D[Zvort,x,x]

Simplify[conc/.ReplaceCStar]

ReplaceAStar=Simplify[Solve[conc==ΛStar/ζ/.ReplaceCStar, A][[1]]]

Geometrical restrictions

By setting the specific values of the geometrical parameters, we obtain the following general restrictions for p=-1

Reduce[Restrictions[-1]/.ReplaceCStar/.ReplaceAStar]

and p=1

Reduce[Restrictions[1]/.ReplaceCStar/.ReplaceAStar]

Thus the set of allowable parameters (p, B, D) can be seen geometrically as

spanBD[Tan2Star_, ΛStar_, p_]=Reduce[FullSimplify[Restrictions[p]/.ReplaceCStar/.ReplaceAStar/.{Γ->1, ζ->1, θ->0}]];
Manipulate[RegionPlot[spanBD[Tan2Star, ΛStar, p], {B, -10, 10}, {D, -10, 10}, FrameLabel->{"B", "D"}, PlotLabel->{Tan2Star, ΛStar, p}], {ΛStar, -1, 1, 0.05}, {Tan2Star, 0.01, 30, 5}, {p, -1, 1, 2}]

Example

An example of set of parameters satisfying the validity conditions is presented below.

GeomValues={Tan2Star->1./3, ΛStar->1, ζStar->1.};
BDPValues={BStar->-5, DStar->0, pStar->-1};
ParamValues={A->-BStar Tan2Star+ΛStar ζStar, B->BStar, C->-BStar Tan2Star+ΛStar ζStar+BStar Tan2Star-DStar Tan2Star, D->DStar, θ->0., p->pStar, ζ->ζStar, Γ->1.}/.BDPValues/.GeomValues;
Restrictions[p]/.ParamValues
If[%, hyp[x_, y_, t_]=Hyperbola/.ParamValues/.{Γ->1};
genICProj=ContourPlot[{hyp[x, y, -1]==0, hyp[x, y, 1]==0}, {x, -2, 2}, {y, -2, 2}, PlotLegends->{"t=-1", "t=1"}, Axes->True, FrameLabel->Automatic, PlotLabel->ParamValues]]

Parametrization of the filaments

Let us now express everything as a function og the geometric parameters of the filaments.

ReplaceTanStar={((Tan[((ϕpStar)/2)])^2)->Tan2Star, ((Cot[((ϕpStar)/2)])^2)->1/Tan2Star};
ReplaceTan2StarInv={Tan2Star->TanStar^2};
ReplacePSquare={p^2->1};

The projections of the filaments onto the z=0 plane and onto the y=0 plane result respectively in

HyperbolaGeom=Simplify[HyperbolaSimp/.ReplaceCStar/.ReplaceAStar/.ReplacePSquare]

HyperbolaGeom=Simplify[HyperbolaSimp/.ReplaceCStar/.ReplaceAStar/.ReplaceTan2StarInv/.ReplacePSquare]

and

ParabolaGeom=FullSimplify[z==Zvort/.Replaceθ0/.ReplaceCStar/.ReplaceAStar/.ReplacePSquare]

ParabolaGeom=FullSimplify[z==Zvort/.Replaceθ0/.ReplaceCStar/.ReplaceAStar/.ReplaceTan2StarInv/.ReplacePSquare]

Following the Supplemental Material, the x, y, and z components of the filaments are parametrized according to

Vort1M={x->-δM/2 1/Tan[((ϕpStar)/2)] Sinh[l], y->δM/2 Cosh[l], z->zM};
Vort2M={x->δM/2 1/Tan[((ϕpStar)/2)] Sinh[l], y->-δM/2 Cosh[l], z->zM};
Vort1P={x->-δP/2 Cosh[l], y->δP/2 Tan[((ϕpStar)/2)] Sinh[l], z->zP};
Vort2P={x->δP/2 Cosh[l], y->-δP/2 Tan[((ϕpStar)/2)] Sinh[l], z->zP};

We can therefore find the functions δM(t) and δP(t) expressed in terms of the initial condition and geometrical parameters.

δMSol=Solve[FullSimplify[HyperbolaGeom/.Vort1M/.ReplaceTanStar/.{l->0},{Restrictions[p]}], δM]//FullSimplify
δPSol=Solve[FullSimplify[HyperbolaGeom/.Vort1P/.ReplaceTanStar/.{l->0},{Restrictions[p]}], δP]//FullSimplify

We can calculate that the ratio AP/AM by taking the square root of the square of the ratio δP/δM, in order to get rid of t and Γ.

Sqrt[((δP/.δPSol[[2]]/.t->Abs[t])/(δM/.δMSol[[2]]/.t->-Abs[t]))^2]//FullSimplify

Finally, we can express the functions zM(l, t) and zP(l, t) in terms of the initial condition, geometrical parameters and x-component of the filaments before and after the reconnection, respectively.

Collect[ParabolaGeom, t]

zMSol=Solve[FullSimplify[ParabolaGeom/.Vort1M/.ReplaceTanStar,{Restrictions[p]}], zM]//FullSimplify
zPSol=Solve[FullSimplify[ParabolaGeom/.Vort1P/.ReplaceTanStar,{Restrictions[p]}], zP]//FullSimplify

Momentum and energy differences during a reconnection

We assume that the filaments are separated by a distance δ: this means setting δM=δP=δ. We want to compute the change in the linear momentum and energy of the filament during the reconnection, occurring within a cylinder of radius R having cylindrical axis parallel equal to the z-axis.

ReplaceAr={Tan[((ϕpStar)/2)]->1/Ar, Cot[((ϕpStar)/2)]->Ar, Tan2Star->1/Ar^2}

R1mAr[l_, δM_, Ar_]={x,y,z}/.Vort1M/.zMSol[[1]]/.ReplaceAr//FullSimplify
R2mAr[l_, δM_, Ar_]={x,y,z}/.Vort2M/.zMSol[[1]]/.ReplaceAr//FullSimplify
R1pAr[l_, δP_, Ar_]={x,y,z}/.Vort1P/.zPSol[[1]]/.ReplaceAr//FullSimplify
R2pAr[l_, δP_, Ar_]={x,y,z}/.Vort2P/.zPSol[[1]]/.ReplaceAr//FullSimplify

The cylindrical region of integration

In terms of the filament parametrization l, the bounds are given by the equations

Lm[Ar_,δM_,R_]=l/.FullSimplify[Solve[R1mAr[l, δM, Ar][[1]]^2+R1mAr[l, δM, Ar][[2]]^2==R^2,l],{δM>0, Ar2>0, R>0}][[4]]/.C[1]->0
Lp[Ar_,δP_,R_]=l/.FullSimplify[Solve[R1pAr[l, δP, Ar][[1]]^2+R1pAr[l, δP, Ar][[2]]^2==R^2,l],{δP>0, Ar2>0, R>0}][[4]]/.C[1]->0

When R is going to infinity, keeping δs and Ar finite, the arguments of the logarithms result in

Limit[((8 (R^2)+(-1+(Ar^2)) ((δM)^2)+2 (Sqrt[(4 (R^2)-((δM)^2)) (4 (R^2)+(Ar^2) ((δM)^2))]))/((1+(Ar^2)) ((δM)^2)))/R^2, R->∞]
Limit[((((δP)^2)+(Ar^2) (8 (R^2)-((δP)^2))+2 (Sqrt[(Ar^2) (4 (R^2)-((δP)^2)) (4 (Ar^2) (R^2)+((δP)^2))]))/((1+(Ar^2)) ((δP)^2)))/R^2, R->∞]

LmlargeR[R_]=Log[R/δM]+Log[(4/(Sqrt[1+(Ar^2)]))]
LplargeR[R_]=Log[R/δP]+Log[(4/(Sqrt[1+(Ar^2)]))]+Log[Ar]

LlargeR[R_]=Log[R]+Log[(4/(Sqrt[1+(Ar^2)] δM))]

Linear momentum difference

The linear momentum of the superfluid can be computed by as line integral of the filaments (see more information in the Supplemental Material). The differential elements of linear momentum before and after the reconnection are given by

dPm[δM_,Ar_]=(Cross[R1mAr[l, δM, Ar],D[R1mAr[l, δM, Ar],l]]+Cross[R2mAr[l, δM, Ar],D[R2mAr[l, δM, Ar],l]])//Simplify
dPp[δP_,Ar_]=(Cross[R1pAr[l, δP, Ar],D[R1pAr[l, δP, Ar],l]]+Cross[R2pAr[l, δP, Ar],D[R2pAr[l, δP, Ar],l]])//Simplify

FullSimplify[-(1+Cot[ϕpStar/2]^2)/Cot[ϕpStar/2]]
Plot[-2 Csc[ϕpStar], {ϕpStar, 0, π}]

Energy difference

Within the cylinder of radius R (using the bounds obtained before on the parametrization of the filaments L), the filament length differential element is given by (the factor 2 is because of the two vortices)

dLm[l_,Ar_,Λ_,δ_, ζ_]=FullSimplify[2Sqrt[D[R1mAr[l, δM, Ar],l].D[R1mAr[l, δM, Ar],l]]]/.{δM->δ, ΛStar->Λ, TanStar->1/Ar}
dLp[l_,Ar_,Λ_,δ_, ζ_]=FullSimplify[2Sqrt[D[R1pAr[l, δP, Ar],l].D[R1pAr[l, δP, Ar],l]]]/.{δP->δ, ΛStar->Λ, TanStar->1/Ar}

Special case Λ=0

In the case of Λ=0, there is an analytical expression for the difference of vortex length before and after the reconnection:

dLm0[l_,Ar_,δ_]=FullSimplify[dLm[l,Ar,0,δ, ζ],{Ar>0,δ>0}]
dLp0[l_,Ar_,δ_]=FullSimplify[dLp[l,Ar,0,δ, ζ],{Ar>0,δ>0}]
ΔL0bis=Integrate[dLp0[l,Ar,δ],{l,-Lpp,Lpp},Assumptions->{Ar>0,δ>0,Lpp>0}]-Integrate[dLm0[l,Ar,δ],{l,-Lmm,Lmm},Assumptions->{Ar>0,Lmm>0}];
ΔL0[Ar_,δ_,R_]=FullSimplify[(ΔL0bis/.{Lpp->Lp[Ar,δ,R],Lmm->Lm[Ar,δ,R]}),{δ>0,Ar>0,R>0}]

Below we plot the length difference is given in the case of Λ=0 for different cylinder radii R/δ. As expected, it is weakly dependent on R.

Plot[{ΔL0[Ar,1, 2],ΔL0[Ar,1, 4], ΔL0[Ar,1, 8], ΔL0[Ar,1, 16]},{Ar,0,4}, AxesLabel->{Ar, "ΔL"}, PlotLegends->{"R=2δ", "R=4δ", "R=8δ", "R=16δ"}]
plotΔL0=Plot[{ΔL0[Ar,1, 2],ΔL0[Ar,1, 4], ΔL0[Ar,1, 16]},{Ar,0,4}, AxesLabel->{MaTeX["A_{\\rm r}", Magnification->1.5], MaTeX["\\Delta \\mathcal{L}", Magnification->1.5]}, PlotLegends->Placed[{MaTeX["R/\\delta=2", Magnification->1.5], MaTeX["R/\\delta=4", Magnification->1.5], MaTeX["R/\\delta=16", Magnification->1.5]}, Bottom], PlotLabel->MaTeX["\\delta_{\\rm lin}=1, \\quad \\Lambda=0", Magnification->1.5], PlotRange->{-2, 2}];

Limit case of large R/δ

In the limit of large R/δ, the above result becomes

ΔL0largeR[Ar_,δ_,R_]=FullSimplify[(ΔL0bis/.{Lpp->LplargeR[R],Lmm->LmlargeR[R]})/.{δM->δ, δP->δ},{δ>0,Ar>0,R>0}]

Below we plot the length difference for large R/δ is given in the case of Λ=0 for different cylinder radii R/δ.

Plot[{ΔL0largeR[Ar,1,10], ΔL0largeR[Ar,1,100], ΔL0largeR[Ar,1,1000]},{Ar,0,4}, AxesLabel->{Ar, "ΔL0largeR"}, PlotLegends->{"R=10δ", "R=100δ", "R=1000δ"}]
plotΔL0LargeR=Plot[{ΔL0largeR[Ar,1,10], ΔL0largeR[Ar,1,100], ΔL0largeR[Ar,1,1000]},{Ar,0,4}, AxesLabel->{MaTeX["A_{\\rm r}", Magnification->1.5], MaTeX["\\Delta \\mathcal{L}", Magnification->1.5]}, PlotLegends->Placed[{MaTeX["R/\\delta=10^1", Magnification->1.5], MaTeX["R/\\delta=10^2", Magnification->1.5], MaTeX["R/\\delta=10^3", Magnification->1.5]}, Bottom], PlotLabel->MaTeX["\\delta_{\\rm lin}=1, \\quad \\Lambda=0", Magnification->1.5], PlotRange->{-2, 2}];

Special case Ar=1

In the case of Ar=1, there is an analytical expression for the difference of vortex length before and after the reconnection:

FullSimplify[dLm[l,1,Λ,δ, ζ],{δ>0, ζ>0}]-FullSimplify[dLp[l, 1,Λ,δ, ζ],{δ>0, ζ>0}]

General case for Λ≠0

In the case of Λ≠0, the integrals have to be performed numerically.

ΔLFull[Ar_,Λ_,δ_,R_, ζ_]:=NIntegrate[dLp[l,Ar,Λ,δ, ζ],{l,-Lp[Ar,δ,R],Lp[Ar,δ,R]}]-NIntegrate[dLm[l,Ar,Λ,δ, ζ],{l,-Lm[Ar,δ,R],Lm[Ar,δ,R]}]

Case, δ=1 and R/δ=2 for different Λ/ζ

Rn=2;
δn=1;
ζn=1;
Plot[{ΔLFull[Ar,0,δn,Rn, ζn], ΔLFull[Ar,1,δn,Rn, ζn], ΔLFull[Ar,2,δn,Rn, ζn], ΔLFull[Ar,3,δn,Rn, ζn]},{Ar,0,4}, AxesLabel->{Ar, "ΔL"}, PlotLegends->{"|Λ|=0", "|Λ|/ζ=1", "|Λ|/ζ=2", "|Λ|/ζ=3"}, PlotLabel->"δ=1, R/δ=2"]
Plot[{ΔLFull[Ar,0,δn,Rn, ζn], ΔLFull[Ar,10,δn,Rn, ζn], ΔLFull[Ar,100,δn,Rn, ζn], ΔLFull[Ar,1000,δn,Rn, ζn]},{Ar,0,4}, AxesLabel->{Ar, "ΔL"}, PlotLegends->{"|Λ|=0", "|Λ|/ζ=10", "|Λ|/ζ=100", "|Λ|/ζ=1000"}, PlotLabel->"δ=1, R/δ=2", PlotRange->All]
plotΔL=Plot[{ΔLFull[Ar,1,δn,Rn, ζn], ΔLFull[Ar,2,δn,Rn, ζn], ΔLFull[Ar,3,δn,Rn, ζn]},{Ar,0,4}, AxesLabel->{MaTeX["A_{\\rm r}", Magnification->1.5], MaTeX["\\Delta \\mathcal{L}", Magnification->1.5]}, PlotLegends->Placed[{MaTeX["|\\Lambda|/\\zeta=1", Magnification->1.5], MaTeX["|\\Lambda|/\\zeta=2", Magnification->1.5], MaTeX["|\\Lambda|/\\zeta=3", Magnification->1.5]}, Bottom], PlotLabel->MaTeX["\\delta_{\\rm lin}=1, \\quad R/\\delta_{\\rm lin}=2", Magnification->1.5]];
plotΔLLargeΛ=Plot[{ΔLFull[Ar,10,δn,Rn, ζn], ΔLFull[Ar,100,δn,Rn, ζn], ΔLFull[Ar,1000,δn,Rn, ζn]},{Ar,0,4}, AxesLabel->{MaTeX["A_{\\rm r}", Magnification->1.5], MaTeX["\\Delta \\mathcal{L}", Magnification->1.5]}, PlotLegends->Placed[{MaTeX["|\\Lambda|/\\zeta=10^1", Magnification->1.5], MaTeX["|\\Lambda|/\\zeta=10^2", Magnification->1.5], MaTeX["|\\Lambda|/\\zeta=10^3", Magnification->1.5]}, Bottom], PlotLabel->MaTeX["\\delta_{\\rm lin}=1, \\quad R/\\delta_{\\rm lin}=2", Magnification->1.5], PlotRange->All];

Case, δ=2 and R/δ=2 for different Λ/ζ

Rn=4;
δn=2;
ζn=1;
Plot[{ΔLFull[Ar,0,δn,Rn, ζn], ΔLFull[Ar,1,δn,Rn, ζn], ΔLFull[Ar,2,δn,Rn, ζn], ΔLFull[Ar,3,δn,Rn, ζn]},{Ar,0,4}, AxesLabel->{Ar, "ΔL"}, PlotLegends->{"|Λ|=0", "|Λ|/ζ=1", "|Λ|/ζ=2", "|Λ|/ζ=3"}, PlotLabel->"R/δ=16"]
Plot[{ΔLFull[Ar,0,δn,Rn, ζn], ΔLFull[Ar,10,δn,Rn, ζn], ΔLFull[Ar,100,δn,Rn, ζn], ΔLFull[Ar,1000,δn,Rn, ζn]},{Ar,0,4}, AxesLabel->{Ar, "ΔL"}, PlotLegends->{"|Λ|=0", "|Λ|/ζ=10", "|Λ|/ζ=100", "|Λ|/ζ=1000"}, PlotLabel->"R/δ=16", PlotRange->All]

Dependence on δ

RRatio=10^3;
ζn=1;
Arn=2;
Plot[{ΔLFull[Arn,0.,δn,RRatio δn, ζn], ΔLFull[Arn,0.5,δn,RRatio δn, ζn], ΔLFull[Arn,1,δn,RRatio δn, ζn], ΔLFull[Arn,2,δn,RRatio δn, ζn]}, {δn, 0.1, 50},  AxesLabel->{"δ", "ΔL"}, PlotLegends->{"|Λ|=0", "|Λ|/ζ=0.5", "|Λ|/ζ=1", "|Λ|/ζ=2"}, PlotLabel->"Ar=2, R/δ=1000"]
plotΔLδΛ=Plot[{ΔLFull[Arn,0.,δn,RRatio δn, ζn], ΔLFull[Arn,0.5,δn,RRatio δn, ζn], ΔLFull[Arn,1,δn,RRatio δn, ζn], ΔLFull[Arn,2,δn,RRatio δn, ζn]}, {δn, 0.1, 50},  AxesLabel->MaTeX@{"\\delta_{\\rm lin}", "\\Delta \\mathcal{L}"}, PlotLegends->Placed[MaTeX@{"|\\Lambda|=0", "|\\Lambda|/\\zeta=0.5", "|\\Lambda|/\\zeta=1", "|\\Lambda|/\\zeta=2"}, Bottom], PlotLabel->MaTeX@"A_{\\rm r}=2, \\quad R/\\delta_{\\rm lin}=10^3"];

RRatio=10^3;
ζn=1;
Λn=1;
Plot[{ΔLFull[0.25,Λn,δn,RRatio δn, ζn], ΔLFull[0.5,Λn,δn,RRatio δn, ζn], ΔLFull[2.,Λn,δn,RRatio δn, ζn], ΔLFull[4,Λn,δn,RRatio δn, ζn]}, {δn, 0.1, 50},  AxesLabel->{"δ", "ΔL"}, PlotLegends->{"Ar=0.25", "Ar=0.5", "Ar=2", "Ar=2.5"}, PlotLabel->"Λ/ζ=1, R/δ=1000"]
plotΔLδAr=Plot[{ΔLFull[0.25,Λn,δn,RRatio δn, ζn], ΔLFull[0.5,Λn,δn,RRatio δn, ζn], ΔLFull[2.,Λn,δn,RRatio δn, ζn], ΔLFull[4,Λn,δn,RRatio δn, ζn]}, {δn, 0.1, 50},  AxesLabel->MaTeX@{"\\delta_{\\rm lin}", "\\Delta \\mathcal{L}"}, PlotLegends->Placed[MaTeX@{"A_{\\rm r}=1/4", "A_{\\rm r}=1/2", "A_{\\rm r}=2", "A_{\\rm r}=4"}, Bottom], PlotLabel->MaTeX@"|\\Lambda|/\\zeta=1, \\quad R/\\delta_{\\rm lin}=10^3"];