Using c Coefficient to Encode First-order PDE Eigenvalue Problem

I want to apply pdeeig to solve an N = 2 PDE eigenvalue problem. For complex-valued u = [ u(1) , u(2) ]:

−( ∂x − i ∂y ) u(2) + a u(1) = λ u(1)

−( ∂x + i ∂y ) u(1) + b u(2) = λ u(2).

Here, a and b are scalar functions of x and y, but I know how to encode them in the form a u using the a coefficient matrix. My question is whether it's possible to encode the first terms in the form − ∇ · ( c ⊗ ∇u ) using the c coefficient tensor.

On that latter page, there is a formula for the components of this vector, but it is slightly ambiguous. Does the notation mean that if c is a function of x and y, then it will be differentiated via chain rule?

If I interpret the formula in this way, it seems like I can use the following trick:

c(1,2,1,2) = − i x − y

c(1,2,2,1) = i x + y

c(2,1,1,2) = i x − y

c(2,1,2,1) = − i x + y,

with all other components of c zero. Since the mixed partials cancel after chain rule expansion, I'm left with exactly the system above. In the function implementation, of course, x and y are the centroids of the triangles in the mesh, and c is represented by a 16-row matrix, as specified.

I would like to know whether this kind of trick is within the specification of the c coefficient, or whether I have interpreted the documentation incorrectly. Moreover, even if this trick works, I'd like to know if perhaps it might lead to numerical issues.

I have made an attempt to encode my problem using this method, with no good results so far. But I can't come up with a quick way to tell if the issues are due to this technical point or more fundamental issues with my physical setup. I'm hoping the community can help me either implicate or eliminate this technicality as an underlying cause.

[Deepak Ramaswamy]

Edwin, let me comment on the first part of the question. I believe, your C coefficient is correct. I ran it through a PDE coefficients extractor I have and I get identical results:

Before commenting on the numerical validity of the approach, I'd like to understand your problem better. Could you let me know the source of your eigenvalue problem? What do these equations represent?

Regards, Deepak

[Edwin]

Thanks for the reply, Deepak. It is helpful to know that this approach at least makes sense to you.

The system considered here is the 2D Dirac equation. There are some details in getting from the full (3+1)-dimensional Dirac equation to a 2D energy eigenvalue equation, which are somewhat analogous to deriving the time-independent Schrodinger equation. The "time-independent" 2D Dirac equation is:

− i ħc ( Sx ∂x + Sy ∂y ) ψ + mc² Sz ψ = E ψ,

where Sx, Sy, and Sz are the (2-by-2) Pauli spin matrices, m is the mass of the particle, and E is the energy eigenvalue corresponding to the energy eigenstate ψ.

Writing out the matrices and using ψ = [ u(1) ; u(2) ]:

− i ħc ( ∂x − i ∂y ) u(2) + mc² u(1) = E u(1)

− i ħc ( ∂x + i ∂y ) u(1) − mc² u(2) = E u(2).

This is essentially the system I originally posted, except that I had simplified matters by removing the prefactor on the first term and renamed a few variables to follow MATLAB conventions. (I also picked a length scale and energy scale to get rid of units.)

Actually, I just now realized that I had also removed the i. However, in my implementation of the c-coefficient, I did multiply by 1i before returning, so I think that's not the issue.

The interesting part is the mc² term, which I want to think of as spatially-dependent. I'm interested in the case where m = 0 whenever x² + y² ≤ R², and some fixed, positive value everywhere else. If we consider the range of energy eigenvalues E < mc², then the particle is forbidden to be outside the circle (since it would have an energy less than its rest mass energy), but at the same time, the particle can freely exist inside the circle (since it is massless).

Physically, I believe this means there should be discrete, normalizable bound states when E < mc², much like in the finite square well for the Schrodinger equation. I expect something oscillatory inside the circle and exponentially decaying outside. Unfortunately, these are not the results that I am getting from my code.

That's a long-winded explanation, and I'm not sure if I was very clear. However, I hope this does provide some insight into the strange form of my c tensor, the functional forms I'm using for the a matrix (representing mass), the ranges of eigenvalues I'm considering, and the kinds of solutions I expect to find (and why).

Let me know if anything I said here tips you off to possible issues. (For example, numerical problems involving the discontinuity in a, an obvious misinterpretation of the physical meaning of the equations, etc.)