Equivariance and Wigner $D$ matrices for Spherical harmonic rotations

Question

I am trying to understand equivariance in machine learning, specially as discussed in the following paper. Claim is that equivariance is when Group symmetry operation, such as rotation, commutes with the hamiltonian (or any other operation of interest for that matter).

In the above paper and extensive literature online claims that Wigner D matrices have the property that their product with Spherical harmoic representation of a position vector is equivalent to rotating the position vector first, and then taking its Spherical harmonic representation. (Sec 4.1.1 in above paper, equation reproduced below)

$$ Y^{(l)}_m (R(g)\hat{r}) = \sum_{m′} D^{(l)}_{mm′}(g) Y^{(l)}_{m′} (\hat{r}) $$ How even in my python example below, I am unable to demonstrate that, what am I not understanding?

import numpy as np
from scipy.special import sph_harm
import spherical
import quaternionic

def cart2sph(x, y, z):
    r = np.sqrt(x**2 + y**2 + z**2)
    theta = np.arctan2(y,x)
    theta = theta if theta < 0 else theta + 2*np.pi
    phi = np.arccos(z/r)
    return r, theta, phi

x=np.random.rand(3)
r,theta,phi = cart2sph(*x)

# Y = 2
Y2 = [
    sph_harm(-2,2,theta,phi),
    sph_harm(-1,2,theta,phi),
    sph_harm( 0,2,theta,phi),
    sph_harm( 1,2,theta,phi),
    sph_harm( 2,2,theta,phi)
]
Y2 = np.array(Y2)

# Additional rotation
rot_theta = 0.2; rot_phi = 0.43
Y2_rotated = [
    sph_harm(-2,2,theta + rot_theta,phi + rot_phi),
    sph_harm(-1,2,theta + rot_theta,phi + rot_phi),
    sph_harm( 0,2,theta + rot_theta,phi + rot_phi),
    sph_harm( 1,2,theta + rot_theta,phi + rot_phi),
    sph_harm( 2,2,theta + rot_theta,phi + rot_phi)
]
Y2_rotated = np.array(Y2_rotated)

# Wigner D matrix
R = quaternionic.array.from_spherical_coordinates(rot_theta, rot_phi)
wigner = spherical.Wigner(2)
D_ = wigner.D(R)
D = D_[-25:].reshape(-5,5)
print(Y2_rotated -  D @ Y2)

Ideally above code shall give answer 0.

Answer 1

Am not sure what is wrong in above code. But I used E3NN library and following code works as expected:

import torch
import e3nn

x1 = torch.tensor([1.0, 0.0, 0.0])

# Generate Real Spherical Harmonics L=2
K1 = 0.5 * torch.sqrt(torch.tensor(15./torch.pi))
K2 = 0.25 * torch.sqrt(torch.tensor(5./torch.pi))
x,y,z = x1; r2 = torch.sum(x1**2)
x,y,z = z,x,y
Yx = torch.tensor([
    K1 * x*y/r2, 
    K1 * y*z/r2, 
    K2 * (3*z**2 - r2)/r2, 
    K1 * x*z/r2, 
    K1 * 0.5 * (x**2 - y**2)/r2]
    )

# Random unitary matrix
R,Res = torch.linalg.qr(torch.rand(3,3))

# Rotated coordinates SH
Rx1 = R @ x1
x,y,z = Rx1; r2 = torch.sum(Rx1**2)
x,y,z = z,x,y
YRx = torch.tensor([
    K1 * x*y/r2, 
    K1 * y*z/r2, 
    K2 * (3*z**2 - r2)/r2, 
    K1 * x*z/r2, 
    K1 * 0.5 *(x**2 - y**2)/r2
    ])

# Generate Wigner D matrix (This and e3nn.o3.wigner_D gives same results)
ir = e3nn.o3.Irreps("1x2e")
D = ir.D_from_matrix(R)

# ∑𝑚′𝐷(𝑙)𝑚𝑚′(𝑔)𝑌(𝑙)𝑚′(𝑟̂ )
DYx = D @ Yx

print(DYx - YRx)
```