g_{\mu \nu}+h_{\mu \nu}$, where $h_{\mu \nu}$ represents the infinitesimal perturbations to the static spherically symmetric background. Each perturbation component can be decomposed into spherical harmonics. Notice that odd-parity modes get a factor $(-1)^{l+1}$ under parity transformation $(\theta,\phi)\rightarrow (\pi-\theta,\pi+\phi)$ while even-parity modes acquire a factor $(-1)^l$. In this paper, we will focus on the odd-parity (axial) modes. We have in this case $$\begin{aligned}
h_{tt}&=0\,,\\
h_{ti}&=\sum_{l,m}h_{0,(l,m)}(t,r)\Bigl[0,\frac{1}{\sin \theta}\sum_{l,m}\partial_\phi Y_l^m,-\sin\theta\sum_{l,m}\partial_\theta Y_l^m\Bigr]\,,\\
h_{ij}&=\sum_{l,m}\Bigl[h_{1,(l,m)}(t,r)(e_1)_{ij}+h_{2,(l,m)}(t,r)(e_2)_{ij}\Bigr]\,,\end{aligned}$$ where the tensor spherical harmonics $(e_{1,2})_{ij}$ are given by $$\begin{aligned}
(e_1)_{ij}=  
\begin{pmatrix}
    0 & \frac{1}{\sin\theta}\partial_\phi Y_l^m & - \sin\theta\partial_\theta Y_l^m \\
    \frac{1}{\sin\theta}\partial_\phi Y_l^m & 0 & 0\\
    - \sin\theta\partial_\theta Y_l^m & 0 & 0
  \end{pmatrix}\end{aligned}$$ and $$\begin{gathered}
(e_2)_{ij}= \left(
  \begin{matrix}
    0 & 0\\ 
    0 & -\frac{1}{\sin\theta}\Bigl[\partial_{\theta\phi}^2-\cot\theta\partial_\phi\Bigr]Y_l^m \\ 
    0 & \frac{1}{2}\Bigl[-\frac{1}{\sin^2\theta}\partial_\phi^2+\cos\theta\partial_\theta+\sin\theta\partial_\theta^2\Bigr]Y_l^m \\ 
  \end{matrix}\right.                
\\
  \left.
  \begin{matrix}
    0 \\ 
    \frac{1}{2}\Bigl[-\frac{1}{\sin^2\theta}\partial_\phi^2+\cos\theta\partial_\theta+\sin\theta\partial_\theta^2\Bigr]Y_l^m\\
    \Bigl[\sin\theta\partial_{\theta\phi}^2-\cos\theta\partial_\phi\Bigr]Y_l^m
  \end{matrix}\right)\,.\end{gathered}$$

This can be written in a more compact way by introducing $E_{ab}\equiv \sqrt{\det\gamma}~\epsilon_{ab}$ where $a, b$ $\in \{\theta,\phi\}$, $\gamma_{ab}$ being the metric on the 2-sphere, and $\epsilon_{ab}$ the totally anti-symmetric symbol with $\epsilon_{\theta\phi}=1$. We have then $$\begin{aligned}
h_{ta}&=\sum_{l,m}h_{0,(lm)}(t,r)E_{ab}\partial^bY_l^m(\theta,\phi)\,,\\
h_{ra}&=\sum_{l,m}h_{1,(lm)}(t,r)E_{ab}\partial^bY_l^m(\theta,\phi)\,,\\
h_{ab}&=\frac{1}{2}\sum_{l,m}h_{2,(lm)}\Bigl[E_a^{~c}\nabla_{cb}Y_l^m(\theta,\phi)+E_b^{~c}\nabla_{ca}Y_l^m(\theta,\phi)\Bigr]\,.\end{aligned}$$ Not all perturbations are physical, we can fix some perturbations to zero because of the invariance under arbitrary differentiable coordinate transformations. We consider an infinitesimal coordinate transformation in terms of a small 