l=|m| contamination when solving EVP "over S2"

[Nicholas Rui]

Hi all,

I'm currently to solve an EVP "over S2". The context is mostly unimportant, but I'm trying to solve for gravity-mode eigenfunctions under strong magnetic fields.

I am trying to do this problem in general on S2, but, because I need to accommodate NCCs in phi, it seems I can unfortunately not use the spherical harmonic basis. Instead, the main difficulty is trying to solve this problem using some sort of Chebyshev basis in theta and somehow enforcing regularity at the poles. I am running into some pretty unusual (but distinctive) issues which I hope you might have some explanation for.

Four scripts (and diagnostic output) are attached to dense-solve the "same problem", but with varying levels of "abstraction" removed. The test case is to try to recover the m=+1 associated Legendre polynomials.

1. "spharm_test_P_noreduce.py": This tries to directly solve the m=+1 generalized Legendre equation using Chebyshev polynomials as a basis:

(1-u^2)*P'' - 2*u*P' - P/(1-u^2) + lambda*P = 0

Afterwards, upon transforming the output eigenvectors to coefficients in terms of normalized associated Legendre polynomials, I find the following (sensible) behavior:

- The eigenvectors are approximately, but not "exactly", orthogonal to each other.

- The eigenvectors are mostly close to single associated Legendre polynomials, with greater mixing between components for higher eigenvalues.

- The eigenvalues are a few percent off.

I roughly understand this behavior: for odd m, associated Legendre polynomials are not analytic at the endpoints. Regularity conditions assert that they generally scale as ~sqrt(1-u^2), which has an infinite slope at the poles.

2. "spharm_test_P.py": This tries to analytically "factor out" ~sqrt(1-u^2), i.e., solve the problem by replacing P with sqrt(1-u^2)*P everywhere and simplifying the equations. This means solving the following:

(1-u^2)*P'' - 4*u*P' - 2*P + lambda*P = 0

This is the mysterious behavior:

- Every other eigenvector is correct (i.e., a single associated Legendre polynomial), but the others (those with l-1 = even) are contaminated specifically by a (often very large) l=m=1 component.

I'm puzzled, since it seems this equation should be equivalent. I'm not sure if I have accidentally introduced a null direction that Dedalus is trying to find.

For lambda = l*(l+1), WolframAlpha says (expectedly) that the above equation possesses closed-form solutions without using hypergeometric functoins, but they tend not to be analytic at the poles (cf. lambda=2, solved by P'=const./(1-u^2)^2), so I hoped that they wouldn't be in the span of the Chebyshev basis (but maybe that's numerically naive). However, the original generalized Legendre equation also has non-analytic basis functions in its span (Legendre Q functions), so I don't know why this feature would be the issue.

3. "spharm_test_PT.py": This is one layer closer to what I'm actually trying to do. Here, I solve for two fields P, T:

λ*P + s**2*T'' - 4*c*T' - 2*T = 0

-s**2*T + c*T + s**2*P' - c*P = 0

Physically, T is some potential function generating a (physical!) vector field over the sphere. There is definitely a gauge freedom, but on S2 that gauge freedom is just a constant offset which in principle only matters for m=0 (but not necessarily in practice?).

The second equation is actually basically trivial and is solved by P=T. This time:

- Half of the eigenvectors are the same eigenvectors as before, with the same l-1=0 contamination.

- The other half return eigenvalues = infinity, and the eigenvectors have P~0 and T which are (sparse) combinations of a few components.

4. "spharm_test_PTF.py": This is even one layer closer to what I'm trying to do. Here, I solve for three fields P, T, F:

λ*P + s**2*T'' - 4*c*T' - 2*T = 0

-s**2*T' + (1-q)*c*T - 1j*q*s**2*c*F' - 1j*(1-q*c**2)*F + s**2*P' - c*P = 0

-s**2*F' + (1-q)*c*F + 1j*q*s**2*c*T' + 1j*(1-q*c**2)*T - 1j*P = 0

Here, F is something like a stream function for the same vector field, with the same caveats about gauge freedom. Here, q can be taken to be a small number (without which there is a singularity error), but it has a real interpretation related to the rotation rate of the star.

Observed behavior:

- 1/3rd of the eigenvalues/vectors are the same as before (with F=0), 2/3rds have eigenvalues = infinity which look like sparse combinations of T and F

--

Additional context:

My actual equations generalize these further by adding terms (related to magnetism) which couple across different m. I'm ultimately trying to solve those equations (for one set of fields P, T, F for each m included -- not relying on the built-in Fourier basis, and with gauge fixed). My main problem at the moment however boils down to those in the previous examples (l=|m| contamination for modes where l-|m| = even).

Preliminary runs where I don't care about regularity at all and solve the problem in Chebyshev(theta) (not Chebyshev(cos theta) to ensure differentiability near poles) seemed to behave okay because "regularity" was approximately enforced by singular prefactors ~sin theta, but for some modes this wasn't exact and some eigenfunctions would fail to vanish at the poles at the level ~1e-3 when they were supposed to.

I know I'm probably stress-testing Dedalus outside a regime where it may be appropriate, but if you have any idea what's going on please let me know.

Thanks,

Nicholas

CODE OUTPUT:

In [11]: run spharm_test_P_noreduce.py
2026-02-05 14:35:58,337 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 0s, Remaining: 0s, Rate: 2.7e+01/s
Number of infinite eigenvalues: 0 / 30
Eigenvalues: [1.98823972e+00-0.j 5.94378575e+00+0.j 1.18283924e+01-0.j
 1.96451529e+01+0.j 2.92824641e+01-0.j 4.08495240e+01+0.j
 5.40341892e+01+0.j 6.92256764e+01+0.j 8.56984692e+01+0.j
 1.04405843e+02+0.j 1.23942593e+02+0.j 1.46146987e+02+0.j
 1.68741427e+02+0.j 1.94569564e+02+0.j 2.20533868e+02+0.j
 2.50197530e+02+0.j 2.79988458e+02-0.j 3.13655163e+02-0.j
 3.47617644e+02-0.j 3.85360146e+02+0.j 4.23654228e+02+0.j
 4.65482053e+02-0.j 5.08146126e+02+0.j 5.54044527e+02+0.j
 6.01064930e+02+0.j 6.51014438e+02+0.j 7.02363310e+02+0.j
 7.56345956e+02-0.j 6.77451157e+04-0.j 7.72231137e+04+0.j]

In [12]: run spharm_test_P.py
2026-02-05 14:36:11,425 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 0s, Remaining: 0s, Rate: 4.5e+01/s
Number of infinite eigenvalues: 0 / 30
Eigenvalues: [  2.-0.j   6.-0.j  12.-0.j  20.-0.j  30.-0.j  42.-0.j  56.-0.j  72.+0.j
  90.-0.j 110.-0.j 132.-0.j 156.-0.j 182.-0.j 210.-0.j 240.-0.j 272.+0.j
 306.-0.j 342.-0.j 380.-0.j 420.-0.j 462.-0.j 506.-0.j 552.-0.j 600.+0.j
 650.-0.j 702.-0.j 756.-0.j 812.-0.j 870.-0.j 930.-0.j]

In [13]: run spharm_test_PT.py
2026-02-05 14:36:35,335 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 0s, Remaining: 0s, Rate: 2.8e+01/s
Number of infinite eigenvalues: 30 / 60
Eigenvalues: [  2.-0.j   6.+0.j  12.+0.j  20.-0.j  30.+0.j  42.+0.j  56.+0.j  72.+0.j
  90.+0.j 110.-0.j 132.+0.j 156.+0.j 182.+0.j 210.-0.j 240.+0.j 272.+0.j
 306.+0.j 342.-0.j 380.+0.j 420.+0.j 462.+0.j 506.+0.j 552.+0.j 600.+0.j
 650.+0.j 702.+0.j 756.+0.j 812.+0.j 870.+0.j 930.+0.j  inf+0.j  inf+0.j
  inf+0.j  inf+0.j  inf+0.j  inf+0.j  inf+0.j  inf+0.j  inf+0.j  inf+0.j
  inf+0.j  inf+0.j  inf+0.j  inf+0.j  inf+0.j  inf+0.j  inf+0.j  inf+0.j
  inf+0.j  inf+0.j  inf+0.j  inf+0.j  inf+0.j  inf+0.j  inf+0.j  inf+0.j
  inf+0.j  inf+0.j  inf+0.j  inf+0.j]

In [14]: run spharm_test_PTF.py
2026-02-05 14:36:41,852 subsystems 0/1 INFO :: Building subproblem matrices 1/1 (~100%) Elapsed: 0s, Remaining: 0s, Rate: 1.9e+01/s
Number of infinite eigenvalues: 60 / 90
Eigenvalues: [  2.00010001-0.j   6.00010002+0.j  12.00010005-0.j  20.00010009+0.j
  30.00010014-0.j  42.0001002 +0.j  56.00010027+0.j  72.00010035+0.j
  90.00010044-0.j 110.00010054+0.j 132.00010065+0.j 156.00010077-0.j
 182.0001009 -0.j 210.00010104+0.j 240.00010119-0.j 272.00010135-0.j
 306.00010152-0.j 342.0001017 -0.j 380.00010189+0.j 420.00010209-0.j
 462.0001023 -0.j 506.00010252+0.j 552.00010275+0.j 600.00010299-0.j
 650.00010324-0.j 702.0001035 -0.j 756.00010377-0.j 812.00010405-0.j
 870.00010657-0.j 929.95221349-0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j          inf+0.j          inf+0.j
          inf+0.j          inf+0.j]

[Nicholas Rui]

As an interesting addition to this story, Joel Ong looked at this problem and found that by artificially modifying the 0,0 coefficients of the matrices in the solver object in a specific way,

M[0,0] += .25

L[0,0] += -0.5

and solving the problem with these matrices (see the attached code), the resulting l-1=0 "contamination" in the l-1=even eigenvectors drops to ~1e-4. We're not really sure why this happens.

[Nicholas Rui]

Sorry -- original message reposted below without weird formatting

--

[Nicholas Rui]

Resolved.

After some discussion with Daniel Lecoanet, it turns out that this is because the Chebyshev basis is normalized s.t. the n=0 Chebyshev polynomial has an extra 1/sqrt(2) prefactor. Ostensibly this is because the usually normalized Chebyshev polynomials normalize to pi/2 for n>0 but pi for n=0.