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 = 4 * 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.