Eigenvalues of the D-t-N map for cuboids

The script for the paper “Spectral properties of the Dirichlet-to-Neumann map for the Helmholtz equation”

Denis S. Grebenkov,  Michael Levitin, Iosif Polterovich (c) 2026

SetOptions[EvaluationNotebook[],ShowCellLabel->False]
Needs["MaTeX`"]

(* ========================================= *)
(* a Dirichlet eigenvalue of a unit interval *)
(* note notational convention (Subscript[λ, ℵ,0]^Dir) = -∞  *) 
(* ========================================= *)

DirEvalueUnitInterval[ℵ_ /;MemberQ[{"a","s"},ℵ], m_Integer/;m>=0]:=If[m==0,-Infinity,Pi^2 If[ℵ=="s", (2m-1)^2, (2m)^2]];

(* ================================================ *)
(* a Dirichlet eigenvalue (Subscript[λ, ℵ,m]^Dir)((Subscript[Q, α])) of a cuboid (Subscript[Q, α ])*)
(* see definition on page 15 of the preprint        *)
(* ================================================ *)

DirEvalueCuboid[αs_,ℵs_,ms_]/;(VectorQ[αs,Positive]&&VectorQ[ℵs,MemberQ[{"a","s"},#]&]&&VectorQ[ms,IntegerQ[#]&&#>=0&]&&Length[αs]==Length[ℵs]==Length[ms]&&Length[αs]>0):=Sum[DirEvalueUnitInterval[ℵs[[j]],ms[[j]]]/αs[[j]]^2,{j,1,Length[αs]}];

(* ========================================= *)
(* Function (Subscript[f, ℵ])(λ) for a unit interval,      *)
(* see definition on page 15 of the preprint *)
(* ========================================= *)

fUnitInterval[Λ_, ℵ_ /;MemberQ[{"a","s"},ℵ]]:=
If[ℵ=="s",
If[Λ<0,Sqrt[-Λ]Tanh[Sqrt[-Λ]/2], -Sqrt[Λ]Tan[Sqrt[Λ]/2]],
If[Λ<0, Sqrt[-Λ]Coth[Sqrt[-Λ]/2], Sqrt[Λ]Cot[Sqrt[Λ]/2]]
];

(* ==================================================== *)
(* Function σ ↦ (1/(α^2))(Subscript[f, ℵ]^(-1))(α σ) for a rescaled interval, *)
(* see page 15 of the preprint                          *)
(* ==================================================== *)

fIntervalInverse[σ_,ℵ_/;MemberQ[{"a","s"},ℵ], m_Integer/;m>=0,α_:1]/; α>0:=Module[{Λ},1/α^2 (Λ/.FindRoot[fUnitInterval[Λ,ℵ]==σ α, {Λ, If[m==1,DirEvalueUnitInterval[ℵ,m]-2,(DirEvalueUnitInterval[ℵ,m-1]+DirEvalueUnitInterval[ℵ,m])/2],  DirEvalueUnitInterval[ℵ,m-1],DirEvalueUnitInterval[ℵ,m]}][[1]])
];

(* =========================================================== *)
(* Function (Subscript[g, ℵ,m])(σ) for cuboid (Subscript[Q, α]), given by the RHS of (3.9) *)
(* =========================================================== *)

gCuboid[σ_, αs_,ℵs_,ms_]/;(VectorQ[αs,Positive]&&VectorQ[ℵs,MemberQ[{"a","s"},#]&]&&VectorQ[ms,IntegerQ[#]&&#>0&]&&Length[αs]==Length[ℵs]==Length[ms]&&Length[αs]>0):=Sum[fIntervalInverse[σ,ℵs[[j]],ms[[j]],αs[[j]]], {j,1,Length[αs]}];

(* ====================================================== *)
(* Parametric plot of one eigenvalue curve ((Subscript[g, ℵ,m])(σ),σ) *)
(* for cuboid (Subscript[Q, α]), with σ∈[range[[1]], range[[2]]]       *)
(* ====================================================== *)

σCurveCuboid[αs_,ℵs_,ms_,range: {_,_},opt:OptionsPattern[]]/;(VectorQ[αs,Positive]&&VectorQ[ℵs,MemberQ[{"a","s"},#]&]&&VectorQ[ms,IntegerQ[#]&&#>0&]&&Length[αs]==Length[ℵs]==Length[ms]&&Length[αs]>0):=ParametricPlot[{gCuboid[σ, αs,ℵs,ms],σ},{σ,range[[1]],range[[2]]},{opt}];

(* ======== *)
(* Examples *)
(* ======== *)

PrintDirEvalueCuboid[αs_,ℵs_,ms_]/;(VectorQ[αs,Positive]&&VectorQ[ℵs,MemberQ[{"a","s"},#]&]&&VectorQ[ms,IntegerQ[#]&&#>=0&]&&Length[αs]==Length[ℵs]==Length[ms]&&Length[αs]>0):=Module[{de=Sum[DirEvalueUnitInterval[ℵs[[j]],ms[[j]]]/αs[[j]]^2,{j,1,Length[αs]}]},
Print["α = ",αs, ", ℵ= ", ℵs, ", m = ", ms, ", Dirichlet eigenvalue = ", de];
de];

αsexample={2 Pi,3 Pi ,2/3 Pi};
deass111=PrintDirEvalueCuboid[αsexample, {"a","s", "s"}, {1,1,1}];
σCurveCuboid[αsexample, {"a","s", "s"}, {1,1,1}, {-25,5}, PlotRange->{{-20,5.5},Automatic}, Epilog->{Line[{{deass111,-25},{deass111,5}}]}]

deass323=PrintDirEvalueCuboid[αsexample, {"a","s","s"}, {3,2,3}];
deass212=PrintDirEvalueCuboid[αsexample, {"a","s","s"}, {3,2,3}-{1,1,1}];
σCurveCuboid[αsexample, {"a","s","s"}, {3,2,3}, {-40,40}, PlotRange->{{-1,deass323+1},Automatic},AxesOrigin->{0,0},Epilog->{Line[{{deass323,-40},{deass323,40}}],Line[{{deass212,-40},{deass212,40}}]}]

(* ======== *)
(* Figure 5 *)
(* ======== *)

tx={-0.75, -0.65, -0.55}; tx={tx, MaTeX[tx]}//Transpose;
ty={0.95, 1.00, 1.05}; ty={ty, MaTeX[ty]}//Transpose;
lg1=MaTeX["g^{-1}_{(\\mathrm{s},\\mathrm{a}), (1,2)}(\\Lambda)"];
lg2=MaTeX["g^{-1}_{(\\mathrm{a},\\mathrm{s}), (1,2)}(\\Lambda)"];
αs={Pi, 27/8 Pi};
ℵs={"s","a"}; ms={1,2};
p1=σCurveCuboid[αs,ℵs, ms,{0.8,1.2},PlotStyle->];
ℵs={"a","s"}; ms={1,2};
p2=σCurveCuboid[αs,ℵs, ms,{0.8,1.2},PlotStyle->];
figcrossrect=Legended[Show[{p1,p2}, PlotRange->{{-0.8,-0.51}, {0.9,1.07}},Axes->True,AxesOrigin->{-0.5,0.9},AxesLabel->MaTeX[{"\\Lambda","\\sigma^{(\\Lambda)}"}],
Ticks->{tx, ty}],Placed[LineLegend[clr[[1;;2]],{lg1,lg2}], {Center, Top}]]