Armónicos esféricos con Python

Armónicos esféricos¶

Autor: Eduardo Martín Calleja

Mi interés por los armónicos esféricos viene motivado, tal como contaba en mi post anterior por ser estos objetos matemáticos esenciales para la resolución de un gran número de problemas, y en especial para el análisis de las anisotropías del fondo cósmico de microondas. En este post estudiaré que son y mostraré ejemplos de cálculo de los mismos con la librería SciPy de Python, así como una serie de representaciones gráficas de los mismos.

Este artículo está escrito íntegramente utilizando el Notebook de Ipython.

Importaciones y referencias¶

In [1]:

%matplotlib inline
from __future__ import division
import scipy as sci
import scipy.special as sp
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm, colors

# Generar un cuadro con versiones de las librerías utilizadas en este notebook
#https://github.com/jrjohansson/version_information
%load_ext version_information
%version_information scipy, numpy, matplotlib

Out[1]:

Software	Version
Python	2.7.9 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython	3.1.0
OS	Linux 3.13.0 55 generic x86_64 with debian jessie sid
scipy	0.15.1
numpy	1.9.2
matplotlib	1.4.3
Fri Jun 19 23:45:59 2015 CEST

Referencias:

Capítulo XI. "Great Balls of PDEs" del curso en internet de Evans Harrell y James Herod: Linear Methods of Applied Mathematics

Capítulo XII del libro Introduction to Partial Differential Equations de Peter J. Olver. Una versión previa de este capítulo se puede descargar de: Partial Differential Equations inThree–Dimensional Space

Un rápido resumen de la definición y propiedades de los armónicos esféricos, junto con algunos gráficos en Clem Pryke's Logbook

Derivación teórica de la definición de los armónicos esféricos¶

Separación de variables¶

Nos proponemos buscar las soluciones de la ecuación de Laplace:

$$\nabla^2 u = 0$$

cuando esta se expresa en coordenadas esféricas. Ya vimos en el post anterior que el problema puede ser formulado de la siguiente forma:

Siendo $(\rho, \theta, \phi)$ el radio, colatitud (de $0$ a $\pi$), y longitud (de $0$ a $2 \pi$) respectivamente. Se trata de encontrar las funciones $u(\rho, \theta, \phi)$ que satisfagan:

$$ \nabla^2 u = \frac{\partial^2 u}{\partial \rho^2} + \frac{2}{\rho} \frac{\partial u}{\partial \rho} +\frac{1}{\rho^2} \frac{\partial^2 u}{\partial \theta^2} + \frac{\cos(\theta)}{\sin(\theta) \rho^2} \frac{\partial u}{\partial \theta} + \frac{1}{\rho^2 \sin^2 (\theta)} \frac{\partial^2 u}{\partial \phi^2} = 0 $$

Para su resolución, aplicaremos el método de separación de variables, buscando de momento soluciones que son el producto de dos funciones, una que depende solo de $\rho$ y otra que depende solo de las variables angulares $(\theta, \phi)$:

$$u(\rho, \theta, \Phi) = R(\rho) Q(\theta, \phi)$$

Sustituyendo, haciendo las derivaciones, y dividiendo por $u = RQ$ se llega a:

$$ 0 = \frac{\nabla^2 u}{u} = \frac{R''}{R} + \frac{2}{\rho} \frac{R'}{R} + \frac{1}{\rho^2 Q} \underset{\nabla^2_S Q}{\underbrace{\left[ \frac{\partial^2 Q}{\partial \theta^2} + \frac{\cos(\theta)}{\sin(\theta)} \frac{\partial Q}{\partial \theta} + \frac{1}{\sin^2(\theta)} \frac{\partial^2 Q}{\partial \phi^2} \right]}} $$

Donde la expresión entre corchetes la reconocemos, de acuerdo con este post como la Laplaciana Esférica, es decir, el operador laplaciano sobre una superficie esférica de radio unidad, que por lo tanto solo depende de las coordenadas $(\theta, \phi)$. A continuación, multiplicando por $\rho^2$ y agrupando variables en cada lado de la ecuación, llegamos a:

$$ \frac{1}{R(\rho)} \left[ \rho^2 R''(\rho) + 2 \rho R'(\rho) \right] = - \frac{\nabla^2_S Q(\theta, \phi)}{Q(\theta, \phi)} $$

Como los dos miembros dependen de diferentes variables, ambos tienen que ser iguales a una constante $\mu$:

$$ \rho^2 R''(\rho) + 2 \rho R'(\rho) - \mu R(\rho) = 0 $$$$ \nabla^2_S Q(\theta, \phi) + \mu Q(\theta, \phi) = 0 $$

Centrandonos de momento en la segunda ecuación, se trata de una ecuación de valores propios: $- \nabla^2_S Q = \mu Q$. Donde $\mu$ es un valor propio del operador $-\nabla^2_S$. Esta ecuación se llama ecuación esférica de Helmholtz. Las soluciones que buscamos serán vectores propios del operador laplaciano esférico. A estos vectores propios $Q(\theta, \phi)$ del laplaciano esférico se les llama armónicos esféricos

Resolución de la ecuación de Helmholtz¶

Si expandimos la ecuación de Helmholz, esta se escribirá:

$$\frac{\partial^2 Q}{\partial \theta^2} + \frac{\cos(\theta)}{\sin(\theta)} \frac{\partial Q}{\partial \theta} + \frac{1}{\sin^2(\theta)} \frac{\partial^2 Q}{\partial \phi^2} + \mu Q = 0 $$

Para su estudio, vamos a aplicar de nuevo la separación de variables, buscando soluciones de la forma:

$$Q(\theta, \phi) = \Theta(\theta) \Phi(\phi) $$

Sustituyendo, dividiendo por $Q = \Theta \Phi$ y multiplicando por $\sin^2(\theta)$ obtenemos:

$$ \frac{1}{\Theta} \left[ \sin^2(\theta) \Theta'' + \sin(\theta) \cos(\theta) \Theta' \right] + \mu \sin^2 \theta = - \frac{\Phi''}{\Phi}$$

Ambos términos deberán igualar a una constante, llamémosle $\lambda$, por lo que la ecuación de Helmholtz esférica se puede dividir en dos ecuaciones, que expresaremos:

$$\sin^2(\theta) \Theta'' + \sin(\theta) \cos(\theta) \Theta' + (\mu \sin^2 \theta - \lambda) \Theta = 0$$$$\Phi'' = - \lambda \Phi $$

Fijémosnos primero en la segunda ecuación. Si buscamos las soluciones como funciones complejas (lo cual va a simplificar enormemente las expresiones trigonométricas) llegamos a soluciones de la forma.

$$\Phi(\phi) = A e^{i \sqrt\lambda \phi} = A \left[ \cos(\sqrt \lambda \: \phi) + i \sin(\sqrt \lambda \: \phi) \right] $$

Siendo $A$ una constante compleja.

Ahora bien, $\Phi(\phi)$ tiene que ser una función periódica de periodo $2 \pi$, esto es, $\Phi(\phi) = \Phi(\phi + 2 \pi)$, y el caso es que, al incrementar $\phi$ en $2 \pi$ el argumento de las funciones seno y coseno se incrementa en $2 \pi \sqrt \lambda$, por lo que la única posibilidad de que se cumpla la periodicidad es que $\sqrt \lambda$ sea un valor entero $m$, positivo o negativo. $\lambda = m^2$, y las soluciones de la componente azimutal de la ecuación de Helmholtz serán el producto de una constante compleja por:

$$\Phi^m(\phi) = e^{i m \phi} $$

Ahora sustituiremos el valor de $\lambda$ en la primera ecuación:

$$\sin^2(\theta) \Theta'' + \sin(\theta) \cos(\theta) \Theta' + (\mu \sin^2 \theta - m^2) \Theta = 0$$

Ahora hacemos el cambio de variables:

$$z = \cos(\theta) \quad , \quad \Theta(\theta) = P(\cos(\theta)) = P(z)$$

Y de esta manera desaparecen las funciones trigonométricas de la ecuación, quedando:

$$(1-z^2) \frac{d^2 P}{d z^2} - 2z \frac{dP}{dz} + \left( \mu - \frac{m^2}{1-z^2} \right) P = 0 $$

Pensemos algunas propiedades de las soluciones de esta ecuación: en primer lugar, $z = \cos(\theta)$, por lo tanto $-1 \leq z \leq 1$, y la función $P(z)$ necesita solo estar definida en el intervalo $[-1, 1]$. Por otro lado, $\Theta(\theta)$ debe estar bien definida en $\theta = 0$ y en $\theta = \pi$ por lo que $P(z)$ debe estar acotada en $z=-1$ y en $z=1$:

$$| P(-1) | < \infty \quad , \quad | P(1) | < \infty$$

Si ahora suponemos $m$ entero no negativo: $m = 0, 1, 2, \dots$, se cumple que (lo admitiremos sin demostración):

• Los únicos valores de $\mu$ para los que la ecuación anterior admite soluciones no triviales (es decir, valores propios) son de la forma $\mu = l(l+1)$ con $l$ entero y $ l \geq |m|$, por lo que resulta la ecuación siguiente, conocida como ecuación diferencial asociada de Legendre:

$$(1-z^2) \frac{d^2 P}{d z^2} - 2z \frac{dP}{dz} + \left( l(l+1) - \frac{m^2}{1-z^2} \right) P = 0 $$

• Las soluciones compatibles con las condiciones anteriores son las llamadas funciones asociadas de Legendre de primer tipo, de grado $l$ y orden $m$, que se suelen designar por $P_l^m(z)$ (también se las conoce como funciones de Ferrers). Existen también las funciones asociadas de segundo tipo $Q_l^m(z)$ que no es preciso considerar ahora.

El caso más sencillo es el correspondiente a $m=0$ en cuyo caso la ecuación queda reducida a la llamada ecuación diferencial de Legendre:

$$(1-z^2) \frac{d^2 P}{d z^2} - 2z \frac{dP}{dz} + l(l+1) P = 0 $$

Y en este caso las correspondientes funciones de Legendre de orden $m=0$ coinciden con los polinomios de Legendre: $P_l^0(z) = P_l(z) \:, \: l=0, 1, \dots$. Volviendo a las funciones $\Theta(\theta)$. Como $\Theta(\theta) = P(\cos \theta)$ obtenemos las funciones siguientes como soluciones de la componente de la solución de la ecuación de Helmholtz que depende solo del ángulo $\theta$:

$$ \mu_l = l(l+1), \quad \Theta_l^m (\theta) = P_l^m ( \cos \theta) \quad \text{para} \: -l \leq m \leq l$$

Armónicos esféricos¶

Estamos ya preparados para obtener las soluciones de la ecuación esférica de Helmholtz

$$ \nabla^2_S Q(\theta, \phi) + \mu Q(\theta, \phi) = 0 $$

Es decir, los valores y vectores propios del operador Lagrangiano sobre la superficie esférica, $-\nabla^2_S Q$. Buscábamos soluciones del tipo:

$$Q(\theta, \phi) = \Theta(\theta) \Phi(\phi) $$

Y combinando los resultados anteriores, obtenemos los armónicos esféricos $Y_l^m(\theta, \phi)$ ($l$ el grado, $m$ el orden), los cuales se suelen multiplicar por una constante de normalización con el fin de conseguir funciones de norma 1 respecto del producto escalar que vamos a definir más adelante:

$$ \text{Valores propios: } \mu_l = l(l+1)$$$$ \text{Funciones propias: } Y_l^m(\theta, \phi) = \sqrt{\frac{2l+1}{4 \pi} \frac{(l-m)!}{(l+m)!}} \; \Theta_l^m(\theta) e^{i m \phi} = P_l^m ( \cos \theta) e^{i m \phi} \quad 0 \leq l < \infty, \; -l \leq m \leq l$$

Ortogonalidad de los armónicos esféricos¶

Esta es la propiedad más interesante de los armónicos esféricos, que va a permitir realizar con ellos desarrollos en serie similares a las series de Fourier. para poder hablar de ortogonalidad necesitamos definir un producto escalar en el espacio $\mathsf L^2(\mathbb S^2)$ de las funciones, en general de valores complejos, de cuadrado integrable sobre la superficie esférica unidad. En la definición se hará uso del elemento diferencial de ángulo sólido en coordenadas esféricas (o si se prefiere, el elemento diferencial de área sobre una superficie esférica de radio 1), el cual viene dado por $d \Omega = \sin (\theta) \; d \theta \; d \phi$.

$$ \langle f,g \rangle = \int_{\phi=0}^{2 \pi} \int_{\theta=0}^\pi f(\theta, \phi) \overline{g(\theta, \phi)} \sin (\theta) \; d \theta \; d \phi $$

Es decir, es la integral del producto de la primera función por la conjugada compleja de la segunda, extendida a la superficie de la esfera.

Por tratarse de funciones propias del operador diferencial del laplaciano esférico, que es auto-adjunto, los armónicos esféricos correspondientes a valores propios distintos (es decir, valores diferentes del grado $l$, serán ortogonales, es decir, su producto escalar será nulo. Pero la propiedad de ortogonalidad va mas allá, en el sentido de que varmónicos esféricos con el mismo grado $l$ pero diferente orden $n$ serán tambien ortogonales. Más aún, debido a la constante de normalización que hemos agregado en su definición, podemos hablar de ortonormalidad:

$$ \langle Y_l^m, Y_{l'}^{m'} \rangle = \delta_{l,l'} \: \delta_{m, m'}$$

Expansión en multipolos esféricos¶

El sistema de armónicos esféricos constituye un conjunto ortonormal completo en $\mathsf L^2(\mathbb S^2)$, por lo que cualquier función escalar de cuadrado integrable sobre la esfera $f(\theta, \phi) : \mathbb S^2 \longrightarrow \mathbb C $ admitirá un desarrollo en serie de la forma:

$$f(\theta, \phi) = \sum_{l=0}^\infty \sum_{m=-l}^l a_{l,m} Y^m_ l(\theta, \phi)$$

La serie converge hacia la función según la norma definida en $\mathsf L^2(\mathbb S^2)$. Los coeficientes se pueden obtener multiplicando escalarmente la serie por cada una de las funciones $Y_l^m$ y utilizando el hecho de que forman un conjunto ortonormal:

$$ a_{l,m} = \langle f, Y_l^m \rangle $$

Estos coeficientes son en general complejos, y se cumplen las siguientes relaciones:

$$ Y_l^{-m} = (-1)^m \overline{Y_l^m} \quad \text{para} \: m>0$$$$ a_{l,-m} = (-1)^m \overline{a_{l, m}} \quad \text{para} \: m>0, \; f \; \text{real}$$

A este tipo de desarrollos se le llama multipolar, indicando el valor de $l$ el orden del multipolo. Concretamente, el término en $a_{0,0}$ corresponde al monopolo, los términos en $ a_{1,-1}, a_{1,0}, a_{1,1}$ corresponden al modo dipolo, los términos para $l=2$ serían el modo cuadripolo, para $l = 3$ tendríamos el modo octupolo, y así sucesivamente.

Los armónicos esféricos en Python Scipy¶

La función que los calcula es en este caso: scipy.special.sph_harm(m,l,phi, theta) . Observese que he cambiado la notación utilizada en la documentación de la función para acomodarla a la que se ha venido usando hasta ahora en estas notas, sustituyendo $n, \theta, \phi$ por $l, \phi, \theta$.

Tanto phi como theta pueden ser arrays.

In [2]:

#Ejemplo de cálculo de la función Y_4^2
l = 4
m = 2
theta, phi = 0.6, 0.75    # Unos valores cualesquiera de ángulos en radianes
Y42 = sp.sph_harm(m, l, phi, theta)
Y42

Out[2]:

(0.028428978328405879+0.40088896207789526j)

Vamos a ver como el valor anterior coincide con la definición, constante de normalización incluida. Por un lado las funciones asociadas de Legendre de primer tipo $P_l^m(z)$ se calculan mediante scipy.special.lpmv(m, l, z)

In [3]:

z = np.cos(theta)
P42 = sp.lpmv(m,l,z)
P42

Out[3]:

9.0104878375269504

A continuación calcularemos la constante de normalización:

In [4]:

f = sci.math.factorial
K_norm = np.sqrt((2*l+1)/(4 * np.pi) * f(l-m)/f(l+m))
K_norm

Out[4]:

0.044603102903819275

Si comparamos el valor obtenido para Y42 con el obtenido aplicando directamente la definición, vemos que coinciden:

In [5]:

K_norm * P42* np.exp(m*phi*1j) == Y42

Out[5]:

True

Ortonormalidad de los armónicos esféricos¶

Vamos a comprobar numéricamente que los armónicos esféricos constituyen una familia ortonormal. Vamos a definir una función que calculará el producto interior $\langle f, g \rangle$ de dos funciones de $\theta, \phi$ que toman valores complejos, definidas sobre $\mathbb S^2$

In [6]:

def dotprod(f,g):
    #Scipy no integra directamente funciones complejas. Hay que descomponerlas en dos integrales:
    #parte real y parte imaginaria
    integrand_r = lambda theta, phi: np.real(f(theta, phi) * np.conj(g(theta, phi)) * np.sin(theta))
    integrand_i = lambda theta, phi: np.imag(f(theta, phi) * np.conj(g(theta, phi)) * np.sin(theta))
    rr = sci.integrate.dblquad(integrand_r, 0, 2 * np.pi,lambda theta: 0, lambda theta: np.pi)[0]
    ri = sci.integrate.dblquad(integrand_i, 0, 2 * np.pi,lambda theta: 0, lambda theta: np.pi)[0]
    if np.allclose(rr,0):
        rr = 0
    if np.allclose(ri,0):
        ri=0
    return rr + ri*1j

In [7]:

#Hagamos una primera prueba con dos funciones definidas arbitrariamente:
f = lambda theta, phi: theta * phi
g = lambda theta, phi: np.exp(theta * 1j)
dotprod(f,g)

Out[7]:

(-15.503138340149906-48.70454551700121j)

In [8]:

# Comprobamos la ortogonalidad de los armónicos esféricos:
# Si (l,m) =! (l',m') el producto interior es cero
Y = lambda l, m, theta, phi: sp.sph_harm(m, l, phi, theta)
f = lambda theta, phi: Y(4,3,theta, phi) 
g = lambda theta, phi: Y(4,2,theta, phi) 
dotprod(f,g)

Out[8]:

0j

In [9]:

# Y si (l,m) = (l',m') el producto interior es uno
f = lambda theta, phi: Y(4,3,theta, phi) 
g = lambda theta, phi: Y(4,3,theta, phi) 
dotprod(f,g)

Out[9]:

(1+0j)

Representación gráfica de los armónicos esféricos¶

A continuación se van a construir varias representaciones gráficas del mismo armónico esférico, con el fin de formarnos una idea lo más aproximada posible de su geometría. Comenzaremos por unas representaciones gráficas habituales en tres dimensiones. Dado que se trata de funciones complejas, representaremos primero su módulo, y a continuación su parte real.

In [10]:

l = 4    #grado del armónico esférico
m = 2    # orden
PHI, THETA = np.mgrid[0:2*np.pi:200j, 0:np.pi:100j] #arrays de variables angulares
R = np.abs(sp.sph_harm(m, l, PHI, THETA)) #Array de valores absolutos de Ymn
#A continuación convertimos a coordenadas cartesianas
# para su representación 3D
X = R * np.sin(THETA) * np.cos(PHI)
Y = R * np.sin(THETA) * np.sin(PHI)
Z = R * np.cos(THETA)

Deseamos colorear la superficie en función del valor de la distancia radial respecto del origen, es decir, de $| Y^m_ l(\theta, \phi)|$. El color lo va a definir por tanto el array R, para lo cual se creará un objeto ScalarMappable con valores RGB siguiendo las indicaciones de esta respuesta en Stackoverflow

In [11]:

N = R/R.max()    # Normalizar R para que los colores del plot cubran todo el rango del colormap.
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'), figsize=(12,10))
im = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=cm.jet(N))
ax.set_title(r'$|Y^2_ 4|$', fontsize=20)
m = cm.ScalarMappable(cmap=cm.jet)
m.set_array(R)    # Asignamos al mappable el array de datos sin normalizar,
                  #para que la escala corresponda a los valores de R
fig.colorbar(m, shrink=0.8);

Y a continuación representaremos solo la parte real de Ymn

In [12]:

l = 4    #grado del armónico esférico
m = 2    # orden
PHI, THETA = np.mgrid[0:2*np.pi:200j, 0:np.pi:100j]
R = sp.sph_harm(m, l, PHI, THETA).real

X = R * np.sin(THETA) * np.cos(PHI)
Y = R * np.sin(THETA) * np.sin(PHI)
Z = R * np.cos(THETA)

#Como R tiene valores negativos, utilizaremos una instancia de Normalize
#según http://stackoverflow.com/questions/25023075/normalizing-colormap-used-by-facecolors-in-matplotlib
norm = colors.Normalize()
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'), figsize=(14,10))
m = cm.ScalarMappable(cmap=cm.jet)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=cm.jet(norm(R)))
ax.set_title('real$(Y^2_ 4)$', fontsize=20)
m.set_array(R)
fig.colorbar(m, shrink=0.8);

Como se ve en la anterior figura, los valores negativos y positivos de $\text{Real}( Y_l^m)$ aparecen diferenciados con claridad, al estar representados los negativos en las gamas de azules, y los positivos en escalas del verde al amarillo, dependiendo los valores de la orientación en el espacio. No obstante, prefiero imaginar los armónicos esféricos como campos escalares (presiones, temperaturas, etc.) definidos sobre la esfera unidad, donde toman valores positivos o negativos que se pueden visualizar como crestas o valles sobre esta superficie. La siguiente gráfica va en esta dirección, ya que se va a sumar el radio de la superficie esférica ($R=1$) al valor de la parte real del armónico esférico. Por lo demás el método de obtención del gráfico es el mismo que antes:

In [13]:

l = 4    #grado del armónico esférico
m = 2    # orden
PHI, THETA = np.mgrid[0:2*np.pi:300j, 0:np.pi:150j]
R = sp.sph_harm(m, l, PHI, THETA).real

s = 1
X = (s*R+1) * np.sin(THETA) * np.cos(PHI)
Y = (s*R+1) * np.sin(THETA) * np.sin(PHI)
Z = (s*R+1) * np.cos(THETA)

norm = colors.Normalize()
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'), figsize=(14,10))
m = cm.ScalarMappable(cmap=cm.jet)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=cm.jet(norm(R)))
ax.set_title('1 + real$(Y^2_ 4)$', fontsize=20)
m.set_array(R)
fig.colorbar(m, shrink=0.8);

No obstante, mi representación preferida es una representación bidimensional con la proyección de Mollweide, la cual tiene la propiedad de preservar la proporción entre las áreas. Considero además que puede ser un ejemplo interesante acerca de como construir un "heatmap" sobre la esfera con una proyección de Mollweide.

In [14]:

# Arrays coordenados para la representación gráfica
x = np.linspace(-np.pi, np.pi, 100)
y = np.linspace(-np.pi/2, np.pi/2, 50)
X, Y = np.meshgrid(x, y)

#Arrays de coordenadas esféricas, derivadas de x, y
#Conversiones necesarias para obtener la representación
#de Mollweide correcta
phi = x.copy()    #copia física
phi[x < 0] = 2 * np.pi + x[x<0]
theta = np.pi/2 - y
PHI, THETA = np.meshgrid(phi, theta)

In [15]:

l = 4    #grado del armónico esférico
m = 2    # orden
SH_SP = sp.sph_harm(m, l, PHI, THETA).real    # solo representaremos la parte real

A continuación haremos una representación en forma de Heatmap sobre la esfera con una proyección del Mollweide. El valor $l-m$ representa las "bandas" horizontales, o líneas de latitud constante, donde el armónico esférico se anula. El orden $m$ representa el número de círculos nulos que pasan por los polos. Al contar estas líneas téngase en cuenta que en una proyección plana de la esfera como la de Mollweide cada arco máximo que pasa por los polos es una semicircunferencia. Pueden verse el resultado con diferentes valores de los parámetros $l$ y $m$ en esta excelente página de visualización interactiva

In [16]:

#Para utilizar símbolos en negritas en Latex en los títulos de matplotlib
#según esta respuesta:
#http://stackoverflow.com/questions/14324477/bold-font-weight-for-latex-axes-label-in-matplotlib
matplotlib.rc('text', usetex=True)
matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]

In [17]:

xlabels = ['$210^\circ$', '$240^\circ$','$270^\circ$','$300^\circ$','$330^\circ$',
           '$0^\circ$', '$30^\circ$', '$60^\circ$', '$90^\circ$','$120^\circ$', '$150^\circ$']

ylabels = ['$165^\circ$', '$150^\circ$', '$135^\circ$', '$120^\circ$', 
           '$105^\circ$', '$90^\circ$', '$75^\circ$', '$60^\circ$',
           '$45^\circ$','$30^\circ$','$15^\circ$']

fig, ax = plt.subplots(subplot_kw=dict(projection='mollweide'), figsize=(10,8))
im = ax.pcolormesh(X, Y , SH_SP)
ax.set_xticklabels(xlabels, fontsize=14)
ax.set_yticklabels(ylabels, fontsize=14)
ax.set_title('real$(Y^2_ 4)$', fontsize=20)
ax.set_xlabel(r'$\boldsymbol \phi$', fontsize=20)
ax.set_ylabel(r'$\boldsymbol{\theta}$', fontsize=20)
ax.grid()
fig.colorbar(im, orientation='horizontal');

Aquí se termina este post, que ya ha quedado muy largo. ¡Hasta la próxima!.