Solving the 2D Schrödinger equation with eigensystem, then verifying orthonormality of eigenfunctions with NIntegrate

Question

I am solving the time-independent 2D Schrödinger equation for an interacting electron and hole in the case of anisotropic electron and hole masses, where the interaction is described a modified form of the Coulomb potential known as the Rytova-Keldysh potential.

The goal is to obtain the eigenenergies and eigenfunctions of the bound electron-hole (exciton) system. The problem is that, while the eigenenergies converge to physically realistic values as the NDEigensystem parameter MaxCellMeasure is reduced, some eigenfunctions are not mutually orthogonal, and furthermore, the result of inner products between eigenfunctions are inconsistent depending on the choice of certain parameters in NIntegrate such as MinRecursion. Let us set the problem up before I dive into what is and is not working.

The Schrödinger equation in this case reads:

$$ -\frac{\hbar^2}{2} \left[ \frac{1}{\mu_x} \nabla^2_x + \frac{1}{\mu_y} \nabla^2_y \right] \psi(x,y) + V(x,y) \psi(x,y) = E \psi(x,y), $$ where $x$ and $y$ refer to the $\textit{relative}$ separation between the electron and hole (that is, I've already performed the coordinate transformation to center-of-mass and relative coordinates of the electron-hole system), with $V(x,y)$ given by: $$ V(x,y) = \frac{\pi k e^2}{(\epsilon_1 + \epsilon_2) \rho_0} \left[ H_0 \left( \frac{\sqrt{x^2 + y^2}}{\rho_0} \right) - Y_0 \left( \frac{\sqrt{x^2 + y^2}}{\rho_0} \right) \right] $$

Previously, I successfully obtained the eigenvalues and eigenfunctions for the exciton in the case of isotropic electron-hole masses, in which case I'd performed a coordinate transformation to 2D polar coordinates, and the numerical problem reduces to a 1D Schrödinger equation. See my previous questions here and here, and also note that I used the coordinate transform $r \to Tan[\xi]$ as provided by user Jens in this thread.

The 2D Schrödinger equation, double coordinate transform $\{x \to Tan[\xi], y \to Tan[\psi]\}$, and problem as implemented in NDEigensystem is:

VKeld[kappa_][d_,chi_][x_,y_] := - (1/(4*chi)) * (StruveH[0,Sqrt[x^2 + y^2 + d^2]/(2*Pi*chi/kappa)] - BesselY[0,Sqrt[x^2 + y^2 + d^2]/(2*Pi*chi/kappa)])

Module[
    {
        nmax = 6, mx = 0.15, my = 1.2, pot = VKeld, kappa = 4.89,
        chi = QuantityMagnitude[Quantity[0.41,"Nanometers"],"BohrRadius"],
        d = 0 (* Want to also solve for spatially separated electrons and holes, but focusing on the direct case for now *),
        eps = 9 * 10^-6 (* This is the smallest value of MaxCellMeasure I can use before my computing cluster crashes due to memory limitations *),
        maxi = 10^12,
        vn = None (* Since I need to perform the reverse coordinate transform back to cartesian coordinates, no point in having Mathematica try to normalize the transformed eigenfunctions *),
        tds,
        tdstrans,
        shift=10,
        stranseps=N[(Pi/2)-$MachineEpsilon],
        evs,efs
    },
    tds = (-(1/2)*((1/mx)*D[f[x, y], {x, 2}, {y, 0}] + (1/my)* D[f[x, y], {x, 0}, {y, 2}]) + pot[kappa][d, chi][x, y]*f[x, y] + shift*f[x, y]);
    tdstrans = Simplify[tds /. {f -> (u[ArcTan[#1], ArcTan[#2]] &)} /. {x -> (Tan[\[Xi]]), y -> (Tan[\[Psi]])}, {(Pi/2) > \[Xi] > -(Pi/2), (Pi/2) > \[Psi] > -(Pi/2)}];
    {evs,efs}=NDEigensystem[
        {
            tdstrans,
            DirichletCondition[u[\[Xi], \[Psi]] == 0,Abs[\[Xi]] == (stranseps) || Abs[\[Psi]] == (stranseps)]
        },
     u[\[Xi], \[Psi]],
     {\[Xi], \[Psi]} \[Element] Rectangle[{-stranseps, -stranseps}, {stranseps, stranseps}],
     nmax,
     Method->{"SpatialDiscretization"->{"FiniteElement",{"MeshOptions"->{"MaxCellMeasure"->eps}}},"Eigensystem"->{"Arnoldi","MaxIterations"->maxi},"VectorNormalization"->vn}
    ];
    {
        evs-shift,
        Head/@efs
    }
]

As mentioned, the above code yields well-converged eigenenergies (the ground state eigenenergy is well-converged for MaxCellMeasure->4*10^-4, but a much finer grid is required for convergence of the excited states).

I will note here that the first few eigenstates calculated by the function above correspond to $x$ and $y$ quantum numbers as: $$ |1> \to (n_x = 0, n_y = 0) \\ |2> \to (n_x = 0, n_y = 1) \\ |3> \to (n_x = 0, n_y = 2) \\ |4> \to (n_x = 0, n_y = 3) \\ |5> \to (n_x = 0, n_y = 4) \\ |6> \to (n_x = 1, n_y = 0) \\ |7> \to (n_x = 0, n_y = 5) $$

This correspondence can be deduced by inspecting the eigenfunctions using Plot3D.

As mentioned previously, the problem I'm having is in the results of the off-diagonal inner products such as $<1|2>$, $<1|3>$, etc. as well as the dipole matrix elements corresponding to one-photon optical transitions for linearly-polarized light, that is, $<1|x|2>$, $<1|y|2>$, etc. Not only do these integrals yield obviously un-physical results in some cases, but some choices of parameters for NIntegrate will yield the correct (or at least, a physically plausible) result, while for other integrals those same parameters will yield an unphysical result.

In general, my process for normalizing the eigenfunctions from NIntegrate and performing subsequent inner products is implemented as:

nintaio[rev_, ref_, j_, k_] := Module[ (* This "all-in-one" function normalizes two eigenfunctions from the list of EF's outputted by NDEigensystem, re-defines the normalized eigenfunctions as separate functions, and then performs three integrals based on the inner products <j|k>, |<j|x|k>|^2, and |<j|y|k>|^2 *)
  {
   refs = {ref[[j]], ref[[k]]}, onctab, oncform, foxtab, foxform, 
   foytab, foyform
  },
  (* Define exactly which EFs will be used to conserve memory *)
  {onctab, foxtab, foytab} = Flatten[
     Table[
      Module[
       {
        (* Initialize some local parameters that we want to overwrite on each iteration of the Table *)
        norms, nef, etr,
        s = 10^ess,
        niopts = {Method -> "GlobalAdaptive", MinRecursion -> minr, 
          MaxRecursion -> maxr, MaxPoints -> maxp}
        },
       (* Calculate the normalization constants for each eigenfunction specified by j and k *)
       norms = Table[
         NIntegrate[
          Conjugate[ refs[[i]][ArcTan[x], ArcTan[y]] ] * 
           refs[[i]][ ArcTan[x], ArcTan[y] ], {x, -s, s}, {y, -s, s}, 
          Evaluate@FilterRules[{niopts}, Options[NIntegrate]]], {i, 
          2}];
       (* Define a new variable for the normalized eigenfunctions *)
             nef = 
        Table[ Function[{x, y}, 
          Evaluate[ 
           refs[[i]][ ArcTan[x], ArcTan[y] ] / 
            Sqrt[ norms[[i]] ] ] ], {i, 2}];
       etr = rev[[k]] - rev[[j]];
       {
        NIntegrate[
         Conjugate[ nef[[2]][x, y] ] * nef[[1]][x, y],
         {x, -s, s}, {y, -s, s}, 
         Evaluate@FilterRules[{niopts}, Options[NIntegrate]]],
        2*mx*etr*(NIntegrate[
            Conjugate[ nef[[2]][x, y] ] * x * nef[[1]][x, y],
            {x, -s, s}, {y, -s, s}, 
            Evaluate@FilterRules[{niopts}, Options[NIntegrate]]]^2),
        2*my*etr*(NIntegrate[
            Conjugate[ nef[[2]][x, y] ] * y * nef[[1]][x, y],
            {x, -s, s}, {y, -s, s}, 
            Evaluate@FilterRules[{niopts}, Options[NIntegrate]]]^2)
        }
       ],
      {ess, {3, 4, 5, 6, 7}},
      {maxp, {Automatic, Infinity}},
      {minr, {Automatic, 10, 100, 1000}},
      {maxr, {10^4,10^6}}
    ],
    {{5}, {1}, {2}, {3}, {4}}]
]

I wrote the above function specifically so that I could more easily diagnose which combinations of Options for NIntegrate yield "good" (e.g. physically realistic) results and which combinations yield "bad" (e.g. physically unrealistic) results.

Let's take a look at the results for $j=1,k=2$:

Orthogonality check, <1|2>

Optical transition matrix element, proportional to |<1|x|2>|^2

Optical transition matrix element, proportional to |<1|y|2>|^2

Let me note here that any result for the optical transition matrix element $>1$ is automatically unphysical (see here for a brief explanation)

We can see that $|1>$ and $|2>$ are orthogonal unless the box size $s$ becomes large, but what I find to be noteworthy is that increasing MinRecursion counterintuitively makes the result worse for $s = 10^6$.

The unexpected effect of $s$ and MinRecursion on the result of NIntegrate can also be seen in the table of results for $<1|y|2>$, which should be non-zero since linearly-polarized light in the $y$-direction should be able to be absorbed by the exciton, changing the quantum number $n_y = 0 \to 1$. **For this integral we also find unphysical results for some material parameters at $s=10^4$ and $s=10^5$, even though those same parameters yielded "good" results for the orthogonality integral.

I could continue to link results for $<1|3>$, etc, but I think the point has been made. The options for NIntegrate yield inconsistent, and sometimes unphysical, results, even after normalizing the raw eigenfunctions myself. Furthermore, some combinations of options yield "good" results for one integral but "bad" results for other integrals. Finally, the effect of the choice of NIntegrate options on the integration result is inconsistent, and at times, counter-intuitive or illogical.

I've experimented with changing the NIntegrate Method to "LocalAdaptive","AdaptiveMonteCarlo", etc., but I haven't found that it makes my results consistently more realistic.

Any ideas or suggestions would be greatly appreciated.

EDIT 4/16:

I realize that my post is quite long and rather open-ended regarding the help I'm looking for, so I'll ask more specific questions here...

1) One possible cause of the inconsistencies with NIntegrate could be the InterpolatingFunctions returned by NDEigensystem. I've read elsewhere that NDEigensystem cannot perform calculations with arbitrary precision, that it is instead locked to $MachinePrecision.

1a) Are there other options to NDEigensystem that affect the InterpolatingFunctions?

1b) In general, are there other options for NDEigensystem that might yield a better overall result besides changing MaxCellMeasure?

1c) Is there a way to perform post-processing on the InterpolatingFunctions that would somehow make their behavior more predictable, or the results of their integration more reliably accurate?

2) Perhaps the double-coordinate-transform between Cartesian and tangential coordinates causes trouble with the InterpolatingFunction. Yesterday I tried working in 2D Polar coordinates so that I'd only have to transform $r \to \tan(\psi)$, $\psi \in [0,\infty); \theta \in [0, 2*\pi]$. However, implementing a PeriodicBoundaryCondition seems to consume a tremendous amount of memory - I have 32 GB RAM and can run the Cartesian-transformed system (whose domain is Rectangle[{-Pi/2,-Pi/2},{Pi/2,Pi/2}], for a total area of $\pi^2$) with a MaxCellMeasure = 9 x 10^-6, but for the case of Polar-transformed coordinates, (with domain Rectangle[{0,0},{Pi/2,2*Pi}], again an area of $\pi^2$), I run out of RAM at onlyMaxCellMeasure = 5 x 10^-4, a severe drop in accuracy. If I remove thePeriodicBoundaryCondition, I can again runNDEigensystemwithMaxCellMeasure = 9 x 10^-6without memory issues. Is there something fundamental aboutPeriodicBoundaryCondition` that causes it to consume significantly more memory? I can provide the code for how I'm implementing the Polar-transformed system, either in this question or in a separate question, if necessary.

3) For the Cartesian-transformed system, the only boundary condition that I'm imposing is the DirichletCondition that the function must go to zero at the boundaries. Since this is a second-order differential equation, I realize that I should have a second boundary condition. Imposing Von Neumann BC, that is, specifying that the first derivative of the function should also go to zero at the boundaries, would seem to make the most sense. I'll note, however, that I did not do this previously when I was only solving the 1D Schodinger equation for isotropic electron/hole masses, yet my results seemed to be fine in that case (e.g. physically accurate, well-converged, etc.). I checked the documentation for NeumannValue but the implementation of that function seemed to be very different than DirichletCondition (and the examples given in the docs only described the implementation in DSolve, NDSolve, etc., which takes PDEs in a different form than NDEigensystem). Can I implement Neumann BC using something like DirichletCondition[D[f[\[Xi],\[Psi]],\[Xi]]==0,\[Xi]==Pi/2] (and a similar statement for the \[Psi] coordinate)? Also, I realize that the derivative will be non-trivial due to the coordinate transformation, but I'm just sketching out the code for now. If that's the correct implementation, can I expect that the eigenfunctions will be more well-behaved?

4) I'll admit that I'm not familiar with the advantages/disadvantages of the specific integration strategies for the "GlobalAdaptive" method of NIntegrate. As I mentioned, though, I tried other Methods such as AdaptiveMonteCarlo and did not find that it made my integration results more consistent. Is it possible, though, that different "GlobalAdaptive" strategies would be particularly well-suited to integrating coordinate-transformed InterpolatingFunctions in a large (that is, large enough to be "infinite") Cartesian domain?