From Discrete to Continuous Spectra

Exploring Spectral Distribution for Schrödinger Operators on Finite and Infinite Intervals

Christopher J. Winfield

We study the distribution of eigenspectra for operators of the form -y''+q(x)y with self-adjoint boundary conditions on both bounded and unbounded interval domains. With integrable potentials q, we explore computational methods for calculating spectral density functions involving cases of discrete and continuous spectra where discrete eigenvalue distributions approach a continuous limit as the domain becomes unbounded. We develop methods from classic texts in ODE analysis and spectral theory in a concrete, visually oriented way as a supplement to introductory literature on spectral analysis. As a main result of this study, we develop a routine for computing eigenvalues as an alternative to NDEigenvalues, resulting in fast approximations to implement in our demonstrations of spectral distribution.

Introduction

We follow methods of the texts by Coddington and Levinson [] and by Titchmarsh [] (both publicly available online via archive.org) in our study of the operator ℒ[y](Overscript[=, def])-y''+q(x)y and the associated problem

ℒ[y]=λ y,

where y'(Overscript[=, def])((dy)/(dx)) on the interval ℐ=(0,∞) with real parameter λ and boundary condition

sin(α) y(0)-cos(α) y'(0)= 0

for fixed α, where 0≤α<π. For continuous q∈(L^1)(ℐ) (the set of absolutely integrable functions on ℐ), we study the spectral function ρ(λ) associated with () and () using two main methods: First, following [], we approximate ρ by step functions associated with related eigenvalue problems on finite intervals (Subscript[I, b])(Overscript[=, def])[0,b] for some sufficiently large positive b; then, we apply asymptotic solution estimates along with an explicit formula for spectral density ((d ρ)/(d λ)) []. For some motivation and clarification of terms, we recall a major application: For certain solutions ψ(x,λ) of () and () and for any f∈(L^2)(ℐ) (the set of square-integrable functions on ℐ), a corresponding solution to () may take the form

f(x)=(∫_0^(+∞))g(λ)ψ(x,λ)dρ(λ)=(∫_0^(+∞))g(λ)ψ(x,λ) ((d ρ(λ))/(d λ))dλ

where

g(λ)=(∫_ℐ)ψ(x,λ) f(x)dx

(in a sense described in Theorem 3.1 of Chapter 9 []); here, g is said to be a spectral transform of f. By way of such spectral transforms, the differential operator ℒ may be represented alternatively in the integral form

(∫_ℐ)ψ(x,λ)ℒ[f(x)]dx=λ g(λ),

where ρ induces a measure by which g∈ (L^2)(ℐ,dρ) (roughly, the set of square-integrable functions when integrated against d ρ) and by which Parseval’s equality holds. Typical examples are the complete set of orthogonal eigenfunctions sin(n π x/b): n=1,2,… for α=(π/2) and the corresponding Fourier sine transform in the limiting case b=+∞ (cf. Chapter 9, Section 1 []).

For a fixed, large finite interval (Subscript[I, b]), we consider the problem (), () along with the boundary condition

cos(β) y(b)-sin(β)y'(b)=0,

(0<β≤π), which together admit an eigensystem with correspondence

(Subscript[λ, k])⟺(Subscript[ψ, k])(x),k=1,2,…,

where the eigenvalues (Subscript[λ, k]) satisfying  (Subscript[Δλ, k]) (Overscript[=, def])(Subscript[λ, k+1])-(Subscript[λ, k])>0 and where the eigenfunctions (Subscript[ψ, k])(x) form a complete basis for (L^2)((Subscript[I, b])). Since the associated spectral function (Subscript[ρ, b])(λ) is a step function with jumps at the various (Subscript[λ, k]), we first estimate these (Subscript[λ, k]) by way of a related equation arising from Prüfer (phase-space) variables and compute the corresponding jumps (Subscript[Δρ, b])((Subscript[λ, k]))=((||(Subscript[ψ, k])||)^(-2)).

Then, we use interpolation to approximate the continuous spectral function ρ(λ) using data from a case of large b at points (Subscript[λ, k]) and using

(Subscript[((dρ(λ))/(dλ))|, λ=(Subscript[λ, k])])≈((Δ (Subscript[ρ, b])((Subscript[λ, k])))/(Δ (Subscript[λ, k])))=(1/(((||(Subscript[ψ, k])||)^2)((Subscript[λ, k+1])-(Subscript[λ, k])))),

imposing the condition ρ(λ)≡0 for all λ⩽0.

We compare our results with those of a well-known formula [] appropriate to our case on ℐ, which we outline as follows: For fixed λ>0, let ψ(x,λ) be the solution to () with boundary values

ψ(0,λ)=cos(α); ψ'(0,λ)=sin(α),

for which the asymptotic formula

ψ(x,λ)=A(λ) cos ((Sqrt[λ])x) +B(λ) sin ((Sqrt[λ])x)+o(1)

holds as x→∞. Then we have

((dρ(λ))/(dλ))=(1/(π(Sqrt[λ])((A^2)(λ)+(B^2)(λ))))

from Section 3.5 [].

Finally, in the last section, we apply the above techniques to extend our study to operators on large domains [-b, b] and on ℝ, where spectral matrices take the place of spectral functions as a matrix analog of spectral transforms on these types of intervals (cf. equation (5.5) []). The techniques are described in detail below, but it is of particular interest that our computations uncover an interesting pattern in a discrete-spectrum case, as we are forced to reformulate our approach according to certain eigen-subspaces involved: our desired spectral approximations are resolved by way of an averaging procedure in forming Riemann sums.

Various sections of Chapters 7–9 [] (see also [] and related articles) present useful introductory discussion applied to material presented in this article; yet, with our focus on equations ()–(), one may proceed given basic understanding of Riemann–Stieltjes integration along with knowledge of ordinary differential equations and linear algebra, commensurate with (say) the use of Eigensystem and NDSolve.

An Eigenvalue Estimator

We compute eigenvalues by first computing solutions θ(x,λ) on (Subscript[I, b])×ℝ to the following, arising from Prüfer variables (equation 2.4, Chapter 8 []):

w'=(cos^2)w+(λ-q(x))(sin^2)w; w(0)=arctan(((y(0))/(y'(0))))=(π/2)-α.

Here, tan(θ)=(y/(y')), where y is a nontrivial solution to (), () and () and θ satisfies

(Subscript[θ(b,l)|, l=(Subscript[λ, k])])=β+(k-1) π

for positive integers k. We interpolate to approximate such solutions as an efficient means to invert () in the variable l. And we use the following function on () throughout this article.

Pruefer[initial_,a_,b_,q_]:=ParametricNDSolveValue[
{
y'[x]==Cos[y[x]]^2+(l-q) Sin[y[x]]^2,
y[a]==initial
},y[b],{x,a,b},{l}]

Consider an example with a=0, b=10, and potential q(x)=(e^(-x)) for parameter l with 0≤l≤L, L=5, in the case α=(π/2), β=π.

α1=Pi/2;a1=0;b1=10;q1=Exp[-x];solution1=Pruefer[ArcTan[Sin[α1],Cos[α1]],a1,b1,q1];

We create an interpolation approximation for eigenvalues (Subscript[λ, k]).

β1=Pi;
L1=5;
θApproximation1=Interpolation[Table[{l,solution1[l]},{l,0,L1,.001}]];
EigenvaluesApproximation1[k_]:=l/.NSolve[θApproximation1[l]==β1+(k-1) Pi,l][[1]]

It is instructive to graphically demonstrate the theory behind this method. Here, we consider the eigenvalues as those values of λ where the graph of w=θ(b,λ) intersects the various lines w=β+(k-1) π as we use Floor to find (Subscript[k, max]) (or MaxIndex), our maximum index k, depending on L.

MaxIndex1=Floor[(θApproximation1[L1]-β1)/Pi]+1;
Plot[{θApproximation1[l],Table[β1+(k-1)Pi,{k,MaxIndex1}]},{l,0,L1},AxesLabel->{λ,w},PlotLabel->w==θ[b,λ]]

We choose these boundary conditions so that we may compare our results with those of NDEigenvalues applied to the corresponding problem () and () using DirichletCondition.

operator[q_]:=-y''[x]+q y[x];
WolframEigenvalues=NDEigenvalues[{operator[q1],DirichletCondition[y[x]==0,True]},y[x],{x,a1,b1},MaxIndex1]

CoddingtonLevinsonEigenvalues=Table[Quiet@EigenvaluesApproximation1[k],{k,MaxIndex1}]

We now compare and contrast the methods in this case. The percent differences of the corresponding eigenvalues are all less than 0.2%, even within our limits of accuracy.

Table[Abs[WolframEigenvalues[[n]]-CoddingtonLevinsonEigenvalues[[n]]]/(Mean[{WolframEigenvalues[[n]],CoddingtonLevinsonEigenvalues[[n]]}]),{n,MaxIndex1}]

In contrast, our interpolation method allows some direct control of which eigenvalues are to be computed, whereas NDEigenvalues (in the default setting) outputs a list up to 39 values, starting from the first. Moreover, our method admits nonhomogeneous boundary conditions, where NDEigenvalues admits only homogeneous conditions, Dirichlet or Neumann.

Spectral Density: Discrete Approximation

We proceed to build our approximate spectral density function ((dρ)/(dλ)) for the problem () and () on ℐ with the same potential q as above. We compute eigenvalues likewise but now on a larger interval [0,b] for b=150 and with nonhomogeneous boundary conditions, say given by α=((3π)/4), β=(π/4) (albeit ρ does not depend on β).

a2=0;b2=150;α2=3Pi/4;β2=Pi/4;L2=1;q2=q1;

We compute eigenvalues via our interpolation method and compute a minimum (Subscript[k, min]) (or MinIndex) as well as a maximum index so as to admit only positive eigenvalues; ρ is supported on [0,∞) and negative eigenvalues result in dubious approximations by (Subscript[ρ, b])(λ).

solution2=Pruefer[ArcTan[Sin[α2],Cos[α2]],a2,b2,q2];
θApproximation2=Interpolation[Table[{l,solution2[l]},{l,0,L2,.0005}]];λ2[k_]:=Quiet[l/.FindRoot[θApproximation2[l]==β2+(k-1) Pi,{l,0,-1,1}][[1]]];
MaxIndex2=Floor[(θApproximation2[L2]-β2)/Pi]+1;
ApproximateEigenvalues2=Table[λ2[k],{k,MaxIndex2}];
MinIndex2=Module[{k},k=1;While[ApproximateEigenvalues2[[k]]<=0,k++];k];

We now compute the values ((Subscript[r, k]^2)/(Subscript[Δλ, k])).

ψSolution2=ParametricNDSolve[{operator[q2]==l y[x],y[0]==Cos[α2],y'[0]==Sin[α2]},y,{x,0,b2},{l}];
Δeigenvalues2[k_]:=ApproximateEigenvalues2[[k+1]]-ApproximateEigenvalues2[[k]];
rSquared2[k_]:=1/(NIntegrate[(y[ApproximateEigenvalues2[[k]]][x]/.ψSolution2)^2,{x,0,b2}]);
DensityFunction2[k_]:=rSquared2[k]/Δeigenvalues2[k];
CoddingtonLevinsonPointPlot2=ListPlot[Table[{ApproximateEigenvalues2[[k]],DensityFunction2[k]},{k,MinIndex2,MaxIndex2-1}]]

Fitting Method

We now apply the method of [] as outlined in equation (). We use FindFit to include data from an interval near the endpoint x=b that includes at least one half-period of the period of the fitting functions sin((Sqrt[λ]) x) and cos((Sqrt[λ]) x).

FitTable2[k_]:=Table[{x,y[ApproximateEigenvalues2[[k]]][x]/.ψSolution2},{x,b2-Min[2 Pi/(Sqrt[ApproximateEigenvalues2[[k]]]),b2/2],b2,Min[2 Pi/(Sqrt[ApproximateEigenvalues2[[k]]]),b2/2]/10}];
coefficient[k_]:=FindFit[FitTable2[k],AA Cos[Sqrt[ApproximateEigenvalues2[[k]]] x]+BB Sin[Sqrt[ApproximateEigenvalues2[[k]]] x],{AA,BB},x];
TInterpolation=Interpolation[Table[{ApproximateEigenvalues2[[k]],1/(Pi Sqrt[ApproximateEigenvalues2[[k]]] ((AA/.coefficient[k])^2+(BB/.coefficient[k])^2))},{k,MinIndex2,MaxIndex2}]];
TPlot=Plot[TInterpolation[x],{x,ApproximateEigenvalues2[[MinIndex2]],ApproximateEigenvalues2[[MaxIndex2]]},PlotRange->Full,PlotStyle->Green];

The function ApproximateEigenvalues may return non-numerical results among the first few, in which case we recommend that either b or L be readjusted or that MinIndex be set large enough to disregard such results.

We now compare our results of the discrete and continuous (asymptotic fit) spectral density approximations.

Show[CoddingtonLevinsonPointPlot2,TPlot,AxesLabel->{λ,((Row[{d,ρ}])/(Row[{d,λ}]))},PlotLabel->"Spectral Density Approximations",Epilog->Inset[Framed[SwatchLegend[{Blue,Green},{"Discrete case","Asymptotic fit"}]],Scaled[{.4,.6}]]]

We compare the results by plotting percent differences, all being less than 0.1%.

ListPlot[Table[{ApproximateEigenvalues2[[k]],Abs[TInterpolation[ApproximateEigenvalues2[[k]]]-DensityFunction2[k]]/Mean[{TInterpolation[ApproximateEigenvalues2[[k]]],DensityFunction2[k]}]},{k,MinIndex2,MaxIndex2-1}],AxesLabel->{λ,"%"},PlotLabel->"Percent Differences in Approximations"]

Check with Exact Calculation

We chose q as above because, in part, the solutions can be computed in terms of well-known BesselI (modified Bessel) functions. Replacing λ by l=λ+i ϵ, for λ,ϵ>0, the solutions are linear combinations of

BesselI[2 i (Sqrt[l]),2 (Sqrt[(e^(-x))])],
BesselI[-2 i (Sqrt[l]),2(Sqrt[(e^(-x))])].

From asymptotic estimates (cf. equation 9.6.7 []), we see that the former is dominant and the latter is recessive as x->+∞ when Im l >0. Then, from Chapter 9 [1], equation 2.13 and Theorem 3.1, we obtain the density function by computing

((d ρ(λ))/(d λ))=((-1)/π)(Underscript[lim , ϵ→(SuperPlus[0])])[(lim_(b→+∞))φ(b,λ+i ϵ)/ψ(b,λ+i ϵ)] (Overscript[=, def]) ((-1)/π)(lim_(ϵ→(SuperPlus[0])))m( λ+i ϵ),

where ψ(x,l) is a solution as above and φ(x,l) is a solution with boundary values φ(0,l)=sin(α), φ'(0,l)=-cos(α). (Here, m is commonly known as the Titchmarsh–Weyl m-function.) In the following code, we produce the density function in exact form by replacing functions from (), the dominant by 1 and the recessive by 0, to compute the inside limit and thereafter simply allowing l=λ to be real.

ψExact2=y[x]/.DSolve[{operator[q2]==l y[x],y[0]==Cos[α2],y'[0]==Sin[α2]},y,x][[1]];
ϕExact2=y[x]/.DSolve[{operator[q2]==l y[x],y[0]==Sin[α2],y'[0]==-Cos[α2]},y,x][[1]];
ratio2=ϕExact2/ψExact2

Density2=((-1)/Pi) Im[ratio2/.{BesselI[2 I (Sqrt[l]),2(Sqrt[(E^(-x))])]->1, BesselI[-2 I (Sqrt[l]),2 (Sqrt[(E^(-x))])]->0}]

We likewise compare the exact formula for the continuous spectrum with the discrete results, noting that the exact graph appears to essentially be the same as that obtained by our asymptotic fitting method (not generally expecting the fits to be accurate for small λ!).

ExactPlot2=Plot[Density2,{l,0,L2},PlotRange->Full,PlotStyle->Black];
Show[CoddingtonLevinsonPointPlot2,ExactPlot2,AxesLabel->{λ,((Row[{d,ρ}])/(Row[{d,λ}]))},PlotLabel->"Spectral Density",Epilog->Inset[Framed[SwatchLegend[{Blue,Black},{"Discrete approximation","Continuous, exact"}]],Scaled[{.45,.65}]]]

Extension to Unbounded Domains: A Proof of Concept

For the operator ℒ we now extend our study to large domains (Subscript[ℐ, b])(Overscript[=, def])[-b,b] in the discrete-spectrum case and to the domain ℝ in the continuous-limit case. We choose an odd function potential of the form q(x)=c x (e^(-r|x|)) for positive constants c, r. We focus on the spectral density associated with specific boundary values at x=0 and an associated pair of solutions to (): namely, we consider expansions in the pair (Subscript[ϕ, 1])(x,λ) and (Subscript[ϕ, 2])(x,λ) such that

(Subscript[ϕ, 1])(0,λ)=1,
(Subscript[ϕ, 1]^′)(0,λ)=0,
(Subscript[ϕ, 2])(0,λ)=0,
(Subscript[ϕ, 2]^′)(0,λ)=1.

We apply the above computational methods to the analytical constructs from Chapter 5 [] in both the discrete and continuous cases. First, for the discrete case, we compute spectral matrices associated with self-adjoint boundary-value problems and the pair as in (): We estimate eigenvalues (Subscript[λ, k]):k=1,2,… for an alternative two-point boundary-value problem on (Subscript[ℐ, b]) for (moderately) large b>0 to compute the familiar jumps of the various components (Subscript[ρ, i,j;b])(λ). These components induce measures that appear in the following form of Parseval’s equality for square-integrable functions f on (Subscript[ℐ, b]) (taken in a certain limiting sense):

(Subscript[g, j])(λ)=(∫_(-b)^b)f(t)(Subscript[ϕ, j])(t,λ)dt,
(∫_(-b)^b)((|f(t)|)^2)d t=(∫_(-∞)^∞)((∑_(i,j=1)^2)(Subscript[g, i])(λ)(Subscript[g, j])(λ) d (Subscript[ρ, i,j;b])(λ))d λ

(real-valued case). Second, we compute the various densities as limits as b→+∞ by the formulas

((d (Subscript[ρ, j,k])(λ))/(d λ))=-(1/π)Im (Subscript[M, j,k])(λ),
(Subscript[M, 1,1])=(1/((SuperMinus[m])-(SuperPlus[m]))),
(Subscript[M, 1,,2])=(Subscript[M, 2,1])=(1/2)((SuperMinus[m])+(SuperPlus[m])) (Subscript[M, 1,1]),
(Subscript[M, 2,2])=(SuperMinus[m])(SuperPlus[m])(Subscript[M, 1,1]),

where (SuperPlus[m])(λ) and (SuperMinus[m])(λ) are certain limits of m-functions, related to equation (), but for our ODE problem on domains [0, +∞) and (-∞, 0], respectively. The densities are computed by procedures more elaborate than (), as discussed later. Then, we compare results of the discrete case like in (), approximating

(Subscript[((d(Subscript[ρ, i,j])(λ))/(d λ))|, λ=(Subscript[λ, k])])≈(((Subscript[Δρ, i,j;b])((Subscript[λ, k])))/(Subscript[Δλ, k])).

Discrete Case

After choosing (self-adjoint) boundary conditions (of which the limits (Subscript[ρ, i,j]) happen to be independent)

y(-b)=sin(α),
y'(-b)=cos(α),
y(b)=sin(β),
y'(b)=cos(β),

on an interval (Subscript[ℐ, b]), we estimate eigenvalues and compute coefficients (Subscript[r, 1;k]), (Subscript[r, 2;k]) from the linear combinations

(Subscript[h, k])(x)=(Subscript[r, 1;k])(Subscript[ϕ, 1])(x,(Subscript[λ, k]))+(Subscript[r, 2;k])(Subscript[ϕ, 2])(x,(Subscript[λ, k]))

for the associated orthonormal (complete) set of eigenfunctions (Subscript[h, k]); k=1,2,3,..., whereby

Δ (Subscript[ρ, i,j;b])((Subscript[λ, k]))=(Subscript[r, i;k])·(Subscript[r, j;k])

(real-valued case). Here, the functions (Subscript[h, k])(x) result by normalizing eigenfunctions (Subscript[ψ, k])(x) satisfying () so that we obtain

(Subscript[r, 1;k]) = (((Subscript[ψ, k])(0))/(||(Subscript[ψ, k])||)),
(Subscript[r, 2;k]) = (((Subscript[ψ, k])'(0))/(||(Subscript[ψ, k])||)).

We are ready to demonstrate. Let us choose q(x)=x (e^(-5|x|)), b =50 π and α=(π/2), β=0 (arbitrary). Much of the procedure follows as above, with minor modification, as we include ParametricNDSolveValue to obtain the values (Subscript[ψ, k])(0) and (Subscript[ψ, k])'(0) (the next result may take around three minutes on a laptop).

L3=2;b3=50 Pi;a3=-b3;α3=Pi/2;β3=0;
q3=x Exp[-5Abs[x]];
solution3=Pruefer[ArcTan[Cos[α3],Sin[α3]],a3,b3,q3];
ψsolution3=ParametricNDSolve[{operator[q3]==l y[x],y[a3]==Sin[α3],y'[a3]==Cos[α3]},y,{x,a3,b3},{l}];
ψ3[l_,x_]:=y[l][x]/.ψsolution3;
θApproximation3=Quiet[Interpolation[Table[{l,solution3[l]},{l,0,L3,.001}]]];
MaxIndex3=Floor[(θApproximation3[L3]-β3)/Pi];
λ3[n_]:=Quiet[l/.NSolve[θApproximation3[l]==β3+(n-1) Pi,l]][[1]]
eigenvalues3=Table[λ3[n],{n,MaxIndex3}];
MinIndex3=Module[{k},
k=1;
While[eigenvalues3[[k]]<=0,k++];k
];
NormSquared3[n_]:=NIntegrate[ψ3[eigenvalues3[[n]],x]^2,{x,a3,b3}];
NormFactor3=Table[NormSquared3[n],{n,MinIndex3,MaxIndex3}];
ψAtZero3:=ParametricNDSolveValue[{operator[q3]==l y[x],y[a3]==Sin[α3],y'[a3]==Cos[α3]},y[0],{x,a3,b3},{l}];
ψPrimeAtZero3:=ParametricNDSolveValue[{operator[q3]==l y[x],y[a3]==Sin[α3],y'[a3]==Cos[α3]},y'[0],{x,a3,b3},{l}];
r1[n_]:=Quiet[Evaluate[ψAtZero3[eigenvalues3[[n]]]]/Sqrt[NormFactor3[[n]]]];
r2[n_]:=Quiet[Evaluate[ψPrimeAtZero3[eigenvalues3[[n]]]]/Sqrt[NormFactor3[[n]]]];

We now approximate the density functions by plotting ((Subscript[(OverTilde[λ]), k, ])(Subscript[𝒟, i,j;k])) where

(Subscript[𝒟, i,j;k])(Overscript[=, def])(Subscript[Δρ, i,j;b])((Subscript[λ, k]))/((Subscript[λ, k+2])-(Subscript[λ, k]))+Δ (Subscript[ρ, i,j;b])((Subscript[λ, k+1]))/((Subscript[λ, k+3])-(Subscript[λ, k+1]))

(for certain K≤(Subscript[k, max])-3) as we compute the difference quotients at the various jumps, over even and odd indices separately, and assign the corresponding sums (Subscript[𝒟, i j;k]) to the midpoints (Subscript[(OverTilde[λ]), k ])of corresponding intervals [(Subscript[λ, k]),(Subscript[λ, k+1])].

DifferenceQuotients11=Table[{Mean[{eigenvalues3[[n]],eigenvalues3[[n+1]]}],r1[n]^2/(eigenvalues3[[n+2]]-eigenvalues3[[n]])+r1[n+1]^2/(eigenvalues3[[n+3]]-eigenvalues3[[n+1]])},{n,MinIndex3,MaxIndex3-3,2}];
DifferenceQuotients22=Table[{Mean[{eigenvalues3[[n]],eigenvalues3[[n+1]]}],r2[n]^2/(eigenvalues3[[n+2]]-eigenvalues3[[n]])+r2[n+1]^2/(eigenvalues3[[n+3]]-eigenvalues3[[n+1]])},{n,MinIndex3,MaxIndex3-3,2}];
DifferenceQuotients12=Table[{Mean[{eigenvalues3[[n]],eigenvalues3[[n+1]]}],r1[n] r2[n]/(eigenvalues3[[n+2]]-eigenvalues3[[n]])+r1[n+1] r2[n+1]/(eigenvalues3[[n+3]]-eigenvalues3[[n+1]])},{n,MinIndex3,MaxIndex3-3,2}];
DiscreteCasePlot=ListPlot[{DifferenceQuotients11,DifferenceQuotients22,DifferenceQuotients12},PlotStyle->{Blue,Orange,Green}];

We give the plots below, in comparison with those of the continuous spectra, and give a heuristic argument in the Appendix as to why this approach works.

Continuous Case

First, we apply the asymptotic fitting method using the solutions (Subscript[ϕ, 1]) and (Subscript[ϕ, 2]). Here, we have to compute full complex-valued formulas for the corresponding m-functions (cf. Section 5.7 []) where a slight modification of the derivation of (SuperPlus[m]), via a change of variables and a complex conjugation, results in (SuperMinus[m]) (See Appendix).

ϕ1solution=ParametricNDSolve[{Evaluate[operator[q3]==l y[x]],y[0]==1,y'[0]==0},y,{x,a3,b3},{l}];
ϕ1[l_,x_]:=y[l][x]/.ϕ1solution;
ϕ2solution=ParametricNDSolve[{Evaluate[operator[q3]==l y[x]],y[0]==0,y'[0]==1},y,{x,a3,b3},{l}];
ϕ2[l_,x_]:=y[l][x]/.ϕ2solution;
Fitψ2Plus[l_]:=FindFit[Table[{x,ϕ2[l,x]},{x,b3-2 Pi/Sqrt[l+1],b3,2 Pi/(15 Sqrt[l+1])}],aa Cos[x Sqrt[l]]+bb Sin[x Sqrt[l]],{aa,bb},x]
Fitψ1Plus[l_]:=FindFit[Table[{x,ϕ1[l,x]},{x,b3-2 Pi/Sqrt[l+1],b3,2 Pi/(15 Sqrt[l+1])}],cc Cos[x Sqrt[l]]+dd Sin[x Sqrt[l]],{cc,dd},x]
Fitψ2Minus[l_]:=FindFit[Table[{x,ϕ2[l,x]},{x,a3,a3+2 Pi/Sqrt[l+1],2 Pi/(15 Sqrt[l+1])}],aa Cos[x Sqrt[l]]+bb Sin[x Sqrt[l]],{aa,bb},x]
Fitψ1Minus[l_]:=FindFit[Table[{x,ϕ1[l,x]},{x,a3,a3+2 Pi/Sqrt[l+1],2 Pi/(15 Sqrt[l+1])}],cc Cos[x Sqrt[l]]+dd Sin[x Sqrt[l]],{cc,dd},x]

APlus[l_]:=aa/.Fitψ2Plus[l];
BPlus[l_]:=bb/.Fitψ2Plus[l];
CPlus[l_]:=cc/.Fitψ1Plus[l];
DPlus[l_]:=dd/.Fitψ1Plus[l];

AMinus[l_]:=aa/.Fitψ2Minus[l];BMinus[l_]:=bb/.Fitψ2Minus[l];CMinus[l_]:=cc/.Fitψ1Minus[l];DMinus[l_]:=dd/.Fitψ1Minus[l];

mPlus[l_]:=-(CPlus[l]+I DPlus[l])/(APlus[l]+I BPlus[l]);
mMinus[l_]:=-(CMinus[l]-I DMinus[l])/(AMinus[l]-I BMinus[l]);

dρ11[l_]:=1/(Pi(mMinus[l]-mPlus[l]));
dρ12[l_]:=.5 (mPlus[l]+mMinus[l]) dρ11[l]; 
dρ22[l_]:=mPlus[l] mMinus[l] dρ11[l]; 

AsymptoticFit=Plot[{Im[dρ11[l]],Im[dρ22[l]],Im[dρ12[l]]},{l,eigenvalues3[[MinIndex3]],L3},
PlotStyle->{Blue,Orange,Green}];

We now compare the result of the discrete and asymptotic fitting methods for the elements (Subscript[𝒟, i,j])(Overscript[=, def])((d(Subscript[ρ, i,j]))/(d λ)).

Show[AsymptoticFit,DiscreteCasePlot,Epilog->Inset[
Text@Framed@Grid[{
{"",(Subscript[𝒟, 11]),(Subscript[𝒟, 22]),(Subscript[𝒟, 12])},
{"Discrete",Style["●",Blue,8],Style["●",Orange,8],Style["●",Green,8]},
{"Asymptotic",Style["———",Blue],Style["———",Orange],Style["———",Green]}
},Alignment->{{Left,Center,Center},Automatic}],Scaled[{0.7,0.167}]],PlotLabel->"Spectral Matrix \!\(\*SubscriptBox[\(𝒟\), \(i\\ j\)]\) - Various Methods",ImageSize->{400,300}]

Appendix

We have deferred some discussion on our use of Quiet, comparison of eigenvalue computations, discrete eigenspace decomposition and Weyl m-functions to this section.

First, we have used Quiet to suppress messages warning that some solutions may not be found. From Chapter 8 [], we expect unique solutions since the functions θ(b,·) are strictly increasing. We have also used Quiet to suppress various messages from ParametricNDSolve and other related functions regarding small values of q(x) to be expected with short-range potentials and large domains.

Second, our formulation of (Subscript[𝒟, i,j]) and the midpoints (Subscript[(OverTilde[λ]), k]) as in () arises from a decomposition of the eigenspace by even and odd indices. We motivate this decomposition by an example plot of the values (Subscript[r, 1,k])·(Subscript[r, 2,k]), where the dichotomous behavior is quite pronounced, certainly for large k.

ListPlot[{Table[r1[n] r2[n],{n,1,MaxIndex3,2}],Table[r1[n] r2[n],{n,2,MaxIndex3,2}]},Epilog->Inset[Framed[SwatchLegend[{Blue,Orange},{"k odd","k even"}]],Scaled[{.15,.15}]],PlotLabel->"\!\(TraditionalForm\`The\\ \*SubscriptBox[\(r\), \(1; \\ k\)]·\*SubscriptBox[\(r\), \(2; \\ k\)]\\ Dichotomy\)",AxesLabel->{k,"\!\(\*SubscriptBox[\(r\), \(1; \\ k\)]\) · \!\(\*SubscriptBox[\(r\), \(2; \\ k\)]\)"}]

We are thus inspired to compute the quotients over even and odd indices separately. Then, we consider, say, a relevant expression from Parseval’s equality: for appropriate Fourier coefficients (Subscript[g, i;k]), i=1,2, associated with respective solutions (Subscript[ϕ, i]), we write

(Underscript[(∑^∞), k=1])(Subscript[g, i;k])(Subscript[g, j;k])Δ (Subscript[ρ, i j;b])((Subscript[λ, k]))=(∑_(k odd))((Subscript[g, i;k])(Subscript[g, j;k])Δ (Subscript[ρ, i j;b])((Subscript[λ, k]))+(Subscript[g, i;k+1])(Subscript[g, j;k+1])Δ (Subscript[ρ, i j;b])((Subscript[λ, k+1])))≈(∑_(k odd)) (Subscript[g, i])((Subscript[(OverTilde[λ]), k]))(Subscript[g, j])((Subscript[(OverTilde[λ]), k])) (Subscript[𝒟, i j;k])Δ (Subscript[λ, k]).

We suppose that (Subscript[λ, k+2])-(Subscript[λ, k])≈2Δ (Subscript[λ, k]) and (Subscript[g, i;k])≈(Subscript[g, i;k+1])≈(Subscript[g, i])((Subscript[(OverTilde[λ]), k])) for the corresponding transforms (Subscript[g, i])(λ) in the limit b→∞. Of course, a rigorous argument is beyond the scope of this article.

Finally, we elaborate on the calculations of the m-functions (SuperPlus[m]) and (SuperMinus[m]): Given the asymptotic expressions

(Subscript[ϕ, 1]^±)(x,λ)=(a^±)(λ) cos((Sqrt[λ])x)+(b^±)(λ) sin((Sqrt[λ])x)+o(1),
(Subscript[ϕ, 2]^±)(x,λ)=(c^±)(λ) cos((Sqrt[λ])x)+(d^±)(λ) sin((Sqrt[λ])x)+o(1),

as x→±∞ (resp.), we follow Section 5.7 of [], making changes as needed, with a modification via complex conjugation (l→λ-i ϵ, say) for (SuperMinus[m]) to arrive at

(m^±)(λ)=-(((c^±)(λ)±i (d^±)(λ))/((a^±)(λ)±i (b^±)(λ))) (resp.).

Acknowledgments

The author would like to thank the members of MAST for helpful and motivating discussions concerning preliminary results of this work in particular and Mathematica computing in general.

References

E. A. Coddington and N. Levinson, Theory of Ordinary Differential Equations, New York: McGraw-Hill, 1955. archive.org/details/theoryofordinary00codd.

E. C. Titchmarsh, Eigenfunction Expansions Associated with Second-Order Differential Equations, 2nd ed., London: Oxford University Press, 1962. archive.org/details/eigenfunctionexp0000titc.

E. W. Weisstein. “Operator Spectrum” from MathWorld–A Wolfram Web Resource. mathworld.wolfram.com/OperatorSpectrum.html.

M. Abramowitz and I. A. Stegun, eds., Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, New York: Wiley, 1972.