How to solve PDEs for function $\Psi(z,\bar{z})$ dependent on independent variables $z,\bar{z}$? (Wirtinger derivatives)

Question

For the past few days, I have been struggling to convey to mathematica to solve a PDE that is in terms of the independent variables $(z,\bar{z})$. I know mathematica supports solving PDEs with respect to $(x,y)$, but that conversion would be very very tedious.

The PDE I am looking at is an eigenvalue equation. I'm looking for the plane wave solutions of $$ (-\Delta^{(H)}+\epsilon \hat{Q})\psi(z,\bar{z})= E \psi(z,\bar{z})$$ where $$\Delta^{(H)}= -(1-|z|^2)^2 \cdot \underbrace{4\frac{\partial}{\partial z}\frac{\partial}{\partial \bar{z}}}_{=\Delta}$$ is called the hyperbolic Laplacian, $\Delta = \nabla^2$ refers to the standard Laplacian in $\mathbb{R}^2$ and $\epsilon>0$ is a small parameter.

Denote by $$\frac{\partial}{\partial z}:= \partial_{z}\ ;\ \frac{\partial}{\partial \bar{z}}:= \bar{\partial_{z}}$$

The Operator $\hat{Q}$ is as follows $$\hat{Q}= \partial_{z}^2 \left((1-|z|^2)^3\partial_{z}\right)+\bar{\partial_{z}}^2 \left((1-|z|^2)^3\bar{\partial_{z}}\right)$$

That is the full problem I'm trying to solve but I wanted to start with a test case first. I am currently struggling to have mathematica consider $z,\bar{z}$ independently yet retain the fact that $z,\bar{z}\in\mathbb{C}$ and then solve the following eigenvalue problem.

Consider the simpler eigenvalue equation for which I know the solutions

$$-\Delta^{(H)} \psi=\varepsilon_{K}\psi$$ With $K=ke^{i\beta}\in\mathbb{C}$, The eigenfunctions have been analytically found to be $$\psi(z)=\left(\frac{1-|z|^2}{|1-ze^{-i\beta}|^2}\right)^{\frac{1+ik}{2}}$$ that satisfy the eigenvalue equation $$-\Delta^{(H)} \psi=\varepsilon_{K}\psi$$ with eigenvalue $\varepsilon_{K}=-(1+k^2)$.

I have verified that this is indeed the eigenfunctions however, It seems I am having trouble telling mathematica to solve this equation as a test case.

I have tried

eqn = -(1 - Abs[z]^2)^2*4*
   D[D[u[z, Conjugate[z]], Conjugate[z]], z];
{vals, funs} = 
  NDEigensystem[{eqn}, 
   u[z, Conjugate[z]], z \[Element] Disk[], 10];

but it is of no avail. I wanted to ask if I could get any help on this matter. Any help would be greatly appreciated!

Answer 1

I will only discuss the test case. An external reference is this PDF by Michael Stone. Some formulas for the Poincare disk model are on this Wikipedia page.

Translation to $x,y$ coordinates. Here $z = x+iy$. Then

DeltaH[f_] := -(1-x^2-y^2)^2*Laplacian[f,{x,y}];

One can see that it is invariant under rotation about the origin, but there is actually a bigger symmetry group, coming from the symmetries of the Poincare disk model, see the PDF above.

Eigenfunctions that OP has given. In $x,y$ coordinates these are

psi[k_,beta_] := ((1-x^2-y^2)/(1+x^2+y^2-2*x*Cos[beta]-2*y*Sin[beta]))^((1+I*k)/2);
lambda[k_,beta_] := -(1+k^2);

Let us check that these are eigenfunctions

- DeltaH[psi[k,beta]] - lambda[k,beta]*psi[k,beta] // Simplify
(* gives 0 *)

The operator is formally selfadjoint on some Hilbert space using an $L^2$-type inner product with measure $dx\,dy / (1-x^2-y^2)^2$. The functions psi[k,beta] are not in that Hilbert space, they have infinite norm. Therefore "eigenfunctions" is to be understood in a loose sense here.

Numerical computation of the eigenvalues. I will use this code:

lambdaN[a_,h_]:=First[NDEigensystem[{
   -(1-x^2-y^2)^a*Laplacian[u[x,y],{x,y}],
   DirichletCondition[u[x,y]==0,True]},
   u[x,y],{x,y}\[Element]Disk[],6,
   Method->{"PDEDiscretization"->{"FiniteElement",
       "MeshOptions"->{"MaxCellMeasure"->h}}}]]//Chop[#,10^(-2)]&;

Here h determines a mesh size, the smaller the better. For DeltaH we must set a=2, but we will also consider the standard Laplacian a=0 and the intermediate value a=1 for comparison.

Standard Laplacian a=0:

lambdaN[0,0.01]
lambdaN[0,0.001]
lambdaN[0,0.0001]
(*
{5.78323,14.6827,14.6827,26.3784,26.3791,30.4787}
{5.78319,14.682,14.682,26.3747,26.3747,30.4713}
{5.78319,14.682,14.682,26.3746,26.3746,30.4713}
*)

The result is quite stable under reduction of mesh size. Compare with known exact values:

{BesselJZero[0,1],BesselJZero[1,1],BesselJZero[2,1],BesselJZero[0,2]}^2//N
(* {5.78319,14.682,26.3746,30.4713} *)

Intermediate case a=1:

lambdaN[1,0.01]
lambdaN[1,0.001]
lambdaN[1,0.0001]
(*
{4.00024,8.00574,8.00574,12.0312,12.0397,16.0613}
{4.00001,8.00017,8.00018,12.001,12.0012,16.002}
{4.,8.,8.,12.,12.,16.0001}
*)

We can again compare with exact results: The eigenfunctions include $ (1-r^2) r^{|k|} e^{ik \varphi}$ with integer $k$, whose eigenvalue is $4(1+|k|)$. We have used polar coordinates $x+iy = re^{i\varphi}$.

The case of interest a=2:

lambdaN[2,0.01]
lambdaN[2,0.001]
lambdaN[2,0.0001]
(*
{2.50649,2.50649,2.97186,3.11274,3.53068,3.53068}
{1.5831,1.95856,1.96664,2.19987,2.22956,2.42164}
{1.66045,1.66154,1.79644,1.79762,1.90218,1.90218}
*)

We get some numbers, but they are not stable under changing the mesh size. This is expected: This operator has continuous spectrum. One would not expect NDEigensystem to give a meaningful result here. Perhaps someone familiar with NDEigensystem can explain how these numbers are computed.

Formulas. Since OP has mentioned this, here are formulas for the derivatives with respect to $z$ in terms of the partial derivatives with respect to $x$ and $y$:

Dz[f_]:=1/2*(D[f,x]-I*D[f,y]);
DzConjugate[f_]:=1/2*(D[f,x]+I*D[f,y]);