spherical harmonics still broken in 9.5.beta8

[Jonathan Thornburg]

Sage has long had problems with spherical harmonics, e.g., this thread
from June 2019:
https://groups.google.com/g/sage-support/c/I_d_meMxRbM/m/Esxo5UO2BAAJ

As of 9.5.beta8, spherical harmonics are (still) broken for some arguments,
with the test case noted in that earlier thread still giving the same
(wrong) result:
sage: theta,phi = var('theta,phi')
sage: spherical_harmonic(1,1,theta,phi)
1/4*sqrt(3)*sqrt(2)*sqrt(sin(theta)^2)*e^(I*phi)/sqrt(pi)
The correct result would be
-1/4*sqrt(6)*e^(I*phi)*sin(theta)/sqrt(pi)
(see, e.g., https://en.wikipedia.org/wiki/Table_of_spherical_harmonics).

In that earlier thread, Eric Gourgoulhon noted that the problem was related
to
https://trac.sagemath.org/ticket/25034
Unfortunately, even though that ticket is now marked "closed defect (fixed)",
spherical harmonics are still broken.

Eric suggested a workaround, using the spin-weighted spherical harmonics
from the /kerrgeodesics_gw/ package, which for spin 0 should reduce to
the standard spherical harmonics.

This works if ell and emm are specific integers, but it can't handle the
case where they are generic variables. E.g., this works with Sage's
native spherical harmonics
sage: ell,emm = var('ell,emm')
sage: theta,phi = var('theta,phi')
sage: diff(spherical_harmonic(ell,emm, theta,phi), theta)
but the analogous
sage: from kerrgeodesic_gw import spin_weighted_spherical_harmonic
sage: diff(spin_weighted_spherical_harmonic(0, ell,emm, theta,phi), theta)
fails.

--
-- "Jonathan Thornburg [remove color- to reply]" <jthor...@pink-gmail.com>
on the west coast of Canada, eh?
"There was of course no way of knowing whether you were being watched
at any given moment. How often, or on what system, the Thought Police
plugged in on any individual wire was guesswork. It was even conceivable
that they watched everybody all the time." -- George Orwell, "1984"

[Eric Gourgoulhon]

Hi,

Le mercredi 5 janvier 2022 à 00:36:11 UTC+1, Jonathan Thornburg a écrit :

Sage has long had problems with spherical harmonics, e.g., this thread
from June 2019:
https://groups.google.com/g/sage-support/c/I_d_meMxRbM/m/Esxo5UO2BAAJ

As of 9.5.beta8, spherical harmonics are (still) broken for some arguments,
with the test case noted in that earlier thread still giving the same
(wrong) result:
sage: theta,phi = var('theta,phi')
sage: spherical_harmonic(1,1,theta,phi)
1/4*sqrt(3)*sqrt(2)*sqrt(sin(theta)^2)*e^(I*phi)/sqrt(pi)
The correct result would be
-1/4*sqrt(6)*e^(I*phi)*sin(theta)/sqrt(pi)
(see, e.g., https://en.wikipedia.org/wiki/Table_of_spherical_harmonics).

 

Actually, the difference between the two results is essentially due to a different convention in the Condon-Shortley phase

(cf. https://en.wikipedia.org/wiki/Spherical_harmonics#Condon%E2%80%93Shortley_phase),

which makes Sage's spherical harmonics Y_l^m differ from Wikipedia and Mathematica ones by a factor (-1)^m.

The other difference in the above example is a lack of simplification of sqrt(sin(theta)^2).

I would vote for including the Condon-Shortley phase in Sage's spherical harmonics, since this is standard in quantum mechanics and this would make Sage agree with Wikipedia and Mathematica.

Eric.

[Eric Gourgoulhon]

Le mercredi 5 janvier 2022 à 08:27:56 UTC+1, Eric Gourgoulhon a écrit :

Actually, the difference between the two results is essentially due to a different convention in the Condon-Shortley phase

(cf. https://en.wikipedia.org/wiki/Spherical_harmonics#Condon%E2%80%93Shortley_phase),

which makes Sage's spherical harmonics Y_l^m differ from Wikipedia and Mathematica ones by a factor (-1)^m.

The other difference in the above example is a lack of simplification of sqrt(sin(theta)^2).

I would vote for including the Condon-Shortley phase in Sage's spherical harmonics, since this is standard in quantum mechanics and this would make Sage agree with Wikipedia and Mathematica.

I've opened

https://trac.sagemath.org/ticket/33117

for this.

In doing so,  I've noticed that current Sage's spherical harmonics disagree with SymPy as well. 

I've also found  a very serious bug in the computation of derivatives of spherical harmonics (see the ticket for details). This has not been seen earlier probably because before https://trac.sagemath.org/ticket/25034 (merged in Sage 9.3), spherical harmonics were basically not usable in Sage.

Eric.

[Jonathan Thornburg]

On Wed, Jan 05, 2022 at 12:14:23AM -0800, Eric Gourgoulhon wrote:
> The other difference in the above example is a lack of simplification of
> sqrt(sin(theta)^2).

I think this problem is worse than "just" a lack of simplification:
if sin(theta) < 0 then sqrt(sin(theta)^2 != sin(theta), i.e., the
theta dependence is wrong, not "just" not-fully-simplified.

Should we have another ticket this?

Eric also wrote:
> Actually, the difference between the two results is essentially due to a
> different convention in the Condon-Shortley phase
> (cf.
> https://en.wikipedia.org/wiki/Spherical_harmonics#Condon%E2%80%93Shortley_phase
> ),
> which makes Sage's spherical harmonics Y_l^m differ from Wikipedia and
> Mathematica ones by a factor (-1)^m.

[[...]]

> I would vote for including the Condon-Shortley phase in Sage's spherical
> harmonics, since this is standard in quantum mechanics and this would make
> Sage agree with Wikipedia and Mathematica.

+1 on this.

> I've opened
> https://trac.sagemath.org/ticket/33117
> for this.

[Eric Gourgoulhon]

Le jeudi 6 janvier 2022 à 06:20:28 UTC+1, Jonathan Thornburg a écrit :

I think this problem is worse than "just" a lack of simplification:
if sin(theta) < 0 then sqrt(sin(theta)^2 != sin(theta), i.e., the
theta dependence is wrong, not "just" not-fully-simplified.

Well, the standard spherical coordinate theta (the colatitude on S^2) lies in [0, pi], so that sin(theta) >= 0.

This simplification issue should be dealt either in https://trac.sagemath.org/ticket/33117 (since it seems pretty easy to implement),  or in a follow up ticket.

Eric. 

[Eric Gourgoulhon]

Hi,

The branch of the ticket https://trac.sagemath.org/ticket/33117 has been merged in Sage 9.6.beta5, so in Sage 9.6 spherical harmonics will agree with those of SymPy, SciPy, Mathematica and Wikipedia, and will have correct derivatives. There remains the issue of simplifying some sqrt(sin(theta)^2) terms which appear for odd orders m. This is now https://trac.sagemath.org/ticket/33501.

Eric.

[Jonathan Thornburg]

Hi Eric,

On Mon, Mar 14, 2022 at 12:58:26PM -0700, Eric Gourgoulhon wrote:
> The branch of the ticket https://trac.sagemath.org/ticket/33117 has been
> merged in Sage 9.6.beta5, so in Sage 9.6 spherical harmonics will agree
> with those of SymPy, SciPy, Mathematica and Wikipedia, and will have
> correct derivatives. There remains the issue of simplifying some
> sqrt(sin(theta)^2) terms which appear for odd orders m. This is now
> https://trac.sagemath.org/ticket/33501.

This is a big improvement already -- thanks a lot for all your efforts
on this (and SageManifolds in general)!

All the best, keep safe and COVID-free, -- Jonathan