Problem solving the Weyl equation

I would like to use FEniCS to solve the Weyl equation, an eigenvalue problem


where σ are the Pauli matrices. When solving the Schroedinger equation -∇^2 u + V u = λ u, I can readily assemble the variational form without much issue>

a = (inner(grad(u), grad(v)) + potential*u*v)*dx                                                                          
m = u*v*dx                                                                                  
A = PETScMatrix()                                                                           
M = PETScMatrix()                                                                           
_ = PETScVector()                                                                           
L = Constant(0.)*v*dx                                                                                                                                                               

assemble_system(a, L, bc, A_tensor=A, b_tensor=_)                                                                                      
assemble_system(m, L, A_tensor=M, b_tensor=_)

But I am not sure how to proceed when the operator ∇^2 is replaced with the operator in the Weyl equation: σ_i . ∂_i

Any help is appreciated.

This looks formally equivalent to the advection equation, which has been studied a lot as a model problem for computational fluid dynamics. In particular, it’s an instance of vector-valued advection–diffusion, i.e., the problem considered by Hughes and Mallet with \mathbf{A}_i = \boldsymbol{\sigma}_i and \mathbf{K} = \mathbf{0}. Solving this robustly usually requires some sort of stabilized scheme, and the linked paper develops the theory for the so-called streamline-upwind Petrov–Galerkin (SUPG) method. It’s not that common to solve the corresponding eigenproblem in fluid mechanics, but I Googled “eigenvalue supg” and found a paper on the (related?) Dirac equation, where the author found that SUPG stabilization suppressed spurious eigenvalues.

The construction of the stabilization term from the first linked paper is rather involved, but I think the math simplifies a lot in the case of \mathbf{A}_i being “Pauli matrices”. Googling those, they look like Hermitian square roots of identity, in which case all of the different variable sets defined by Hughes and Mallet are the same, so I believe the formulation just reduces to: Find \lambda and \mathbf{U} such that, for all \mathbf{V},

\int_\Omega \boldsymbol{\sigma}_i\mathbf{U}_{,i}\cdot\mathbf{V}\,d\Omega + \int_{\Omega'}\tau\left(\boldsymbol{\sigma}_i\mathbf{U}_{,i} - \lambda\mathbf{U}\right)\cdot\left(\boldsymbol{\sigma}_j\mathbf{V}_{,j}\right)\,d\Omega' = \lambda\int_{\Omega}\mathbf{U}\cdot\mathbf{V}\,d\Omega\text{ ,}

where summation over repeated indices is implied, \int_{\Omega'} refers to a sum over cell integrals (i.e., what the dx measure in UFL already does, but the mathematical distinction is technically needed for non-smooth function spaces), and \tau is a mesh-dependent stabilization parameter, of the form \tau := h/2, with h the cell diameter, i.e., h = CellDiameter(mesh) in UFL. (In general, \tau would be matrix-valued, following the construction of Hughes and Mallet, but, again, I’ve attempted to simplify based on my (mis)interpretation of the Wikipedia article about the Weyl equation.)

On assembly, one obtains a generalized eigenvalue problem: Find \lambda and \mathbf{x} such that

\mathbf{A}\mathbf{x} = \lambda\mathbf{B}\mathbf{x}\text{ ,}


\mathbf{A} = \text{assemble}\left(\int_\Omega \boldsymbol{\sigma}_i\mathbf{U}_{,i}\cdot\mathbf{V}\,d\Omega + \int_{\Omega'}\tau\left(\boldsymbol{\sigma}_i\mathbf{U}_{,i}\right)\cdot\left(\boldsymbol{\sigma}_j\mathbf{V}_{,j}\right)\,d\Omega'\right)


\mathbf{B} = \text{assemble}\left(\int_{\Omega}\mathbf{U}\cdot\mathbf{V}\,d\Omega + \int_{\Omega'}\tau\mathbf{U}\cdot\left(\boldsymbol{\sigma}_i\mathbf{V}_{,i}\right)\,d\Omega'\right)\text{ .}