Spherical Harmonics: special cases and more examples

[mjungmath]

This ticket is devoted to improve spherical harmonics and make them function properly. That entails treating special cases more efficiently and adding a bunch of doctests for checking correctness.

Depends on #25034
Depends on #31637

CC: @egourgoulhon @tscrim @slel @mkoeppe

Component: misc

Author: Michael Jung

Branch/Commit: public/spher_harmonics_improved_31639 @ e92db1d

Issue created by migration from https://trac.sagemath.org/ticket/31639

[mjungmath]

Dependencies: #25034

[mjungmath]

Author: Michael Jung

[mjungmath]

Branch: public/spher_harmonics_improved_31639

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

4fcc574	Trac #25034: formula for n=m fixed
a4edb24	Trac #25034: abs(m) exceeding n gives zero immediately
1ace9fa	Trac #31639: special cases added + list of first spherical harmonics added

[sagetrac-git]

Commit: 1ace9fa

[mjungmath]

comment:4

Here is one first approach based on the current branch of #25034.

Any chance to get rid of terms like sqrt(sin^2(theta)) in a nice way?

[sagetrac-git]

Changed commit from 1ace9fa to e92db1d

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

d5461d4	Trac #31639: import statements given when needed
e92db1d	Trac #31639: separate special values from general formula

[egourgoulhon]

comment:7

Replying to @mjungmath:

Here is one first approach based on the current branch of #25034.

Thanks!

Any chance to get rid of terms like sqrt(sin^2(theta)) in a nice way?

These ones are really annoying in Sage; this is why I introduced the functions simplify_sqrt_real and simplify_abs_trig, which are used in the default simplification chain for calculus on manifolds. They do the job here:

sage: x, y = var('x y')
sage: s = spherical_harmonic(2, 1, x, y); s
1/4*sqrt(6)*sqrt(5)*sqrt(sin(x)^2)*cos(x)*e^(I*y)/sqrt(pi)
sage: assume(0 <= x, x <= pi)                              
sage: s.simplify_full()  # no effect despite the above assumptions!
1/4*sqrt(6)*sqrt(5)*sqrt(sin(x)^2)*cos(x)*e^(I*y)/sqrt(pi)

sage: from sage.manifolds.utilities import simplify_sqrt_real, simplify_abs_trig
sage: simplify_abs_trig(simplify_sqrt_real(s))
1/4*sqrt(5)*sqrt(3)*sqrt(2)*cos(x)*e^(I*y)*sin(x)/sqrt(pi)

[egourgoulhon]

comment:8

The following results given in the EXAMPLES section:

Y_1^-1(x,y)=-1/2*sqrt(3/2)*e^(-I*y)*sin(x)/sqrt(pi)
Y_1^1(x,y)=1/4*sqrt(3)*sqrt(2)*sqrt(sin(x)^2)*e^(I*y)/sqrt(pi)
Y_2^-1(x,y)=-1/4*sqrt(15)*sqrt(2)*sqrt(sin(x)^2)*cos(x)*e^(-I*y)/sqrt(pi)
Y_2^1(x,y)=1/4*sqrt(6)*sqrt(5)*sqrt(sin(x)^2)*cos(x)*e^(I*y)/sqrt(pi)

have signs opposite to https://en.wikipedia.org/wiki/Table_of_spherical_harmonics and https://mathworld.wolfram.com/SphericalHarmonic.html. This is certainly due to the Condon-Shortley phase (-1)^m. I think we should agree with Wikipedia's convention, since it is standard in mathematical physics. Problably, we should offer the option condon_shortley=True in spherical_harmonic's arguments. For instance, this is done in the package kerrgeodesic_gw:

sage: from kerrgeodesic_gw import spin_weighted_spherical_harmonic
sage: x, y = var('x y')
sage: spin_weighted_spherical_harmonic(0, 1, -1, x, y)  # weight = 0 -> spherical harmonic
0.25*sqrt(6)*e^(-I*y)*sin(x)/sqrt(pi)
sage: spin_weighted_spherical_harmonic(0, 1, 1, x, y)
-1/4*sqrt(6)*e^(I*y)*sin(x)/sqrt(pi)
sage: spin_weighted_spherical_harmonic(0, 1, 1, x, y, condon_shortley=False)
1/4*sqrt(6)*e^(I*y)*sin(x)/sqrt(pi)

By the way, there is not the sqrt(sin(x)^2) issue discussed in comment:7 in the above outputs, maybe because kerrgeodesic_gw uses the algorithm of Golberg (1967) and does not rely explicitely on associated Legendre function; see the function _compute_sw_spherical_harm in line 27 of the file spin_weighted_spherical_harm.py.

PS: kerrgeodesic_gw can be installed in Sage with sage -pip install kerrgeodesic_gw, but if you want to play with it without installing it, you can run this notebook in Binder.

[egourgoulhon]

comment:9

Another motivation for using the Condon-Shortley phase: SymPy uses it:

sage: x, y = var('x y')                             
sage: from sympy import Ynm                   
sage: Ynm(1, 1, x, y).expand(func=True)
-sqrt(6)*exp(I*y)*sin(x)/(4*sqrt(pi))

We should certainly agree with SymPy (in addition to Wikipedia, MathWorld, etc.)

[mjungmath]

comment:10

Thank you for the input. I agree that a flag for the Condon-Shortley phase is the most favorable option.

I'll probably use simplify_abs_trig. Why haven't these simplifications made it to the global namespace though? I think it is definitely worth to add.

Let's wait for #25034 before we wrap this up here.

[mjungmath]

Changed dependencies from #25034 to #25034, #31637