$N$-layer graphene $J^N_j = \sigma_j^N E_j$. Therefore, the layered conductivity fits the framework of the transfer matrix method.

Because the velocity operator $\hat{v}_j$ is written in terms of Pauli matrices $\sigma_x$ or $\sigma_y$ for pseudospin depending on the light polarization direction, the current matrix element becomes $$\begin{aligned}
 |\langle \Psi_{r{\bf p}}^{s}|(-e \hat{v}^\pm_j) |\Psi_{r'{\bf p}}^{s'}\rangle|^2
 =e^2 (N_{r'{\bf p}}^{s'})^2 (N_{r{\bf p}}^{s})^2 
 \sin^2(jk_r) \sin^2(jk_{r'}) 
 \left|\frac{\Pi}{\varepsilon_{r'{\bf p}}^{s'}} \pm \frac{\Pi^*}{\varepsilon_{r{\bf p}}^{s}} \right|^2,\end{aligned}$$ where the $\pm$ signs denote $x$ and $y$ polarization. The dependence of the matrix element on the polarization is eventually lost for the dynamical conductivity, by the integral about polar angle of ${\bf p}$. Putting above into Eq. (\[eq:layeredsigma\]) gives $$\begin{aligned}
 \sigma^N_j(\omega) &= g_s g_v \frac{e^2}{i\hbar} \sum_{s,s'} \frac{1}{2\pi} \int pdp \sum_{r,r'}
 \frac{f(\varepsilon_{r'{\bf p}}^{s'})-f(\varepsilon_{r{\bf p}}^{s})}{(\varepsilon_{r'{\bf p}}^{s'}-\varepsilon_{r{\bf p}}^{s})(\varepsilon_{r'{\bf p}}^{s'}-\varepsilon_{r{\bf p}}^{s}+\hbar \omega + i\epsilon)} \nn \\
 & \times \frac{1}{sE_{r}} \frac{\sin^2(jk_r)}{N+1} \times
 \frac{1}{s'E_{r'}} \frac{\sin^2(jk_{r'})}{N+1} \times
 p^2 \left(\frac{\varepsilon_{r{\bf p}}^{s}}{\varepsilon_{r'{\bf p}}^{s'}} + \frac{\varepsilon_{r'{\bf p}}^{s'}}{\varepsilon_{r{\bf p}}^{s}} \right).\end{aligned}$$ We evaluate this for charge neutral condition $E_F=0$ and at zero temperature $T=0$. Then, the Fermi distribution functions become the step functions. As a result, since $\varepsilon_{r{\bf p}}^{+} \ge 0$ and $\varepsilon_{r{\bf p}}^{-} \le 0$, only interband transitions contribute to the dynamical conductivity; $$\begin{aligned}
 \sigma^N_j(\omega)
 &= -g_s g_v \frac{e^2}{i\hbar} \frac{1}{2\pi} \int pdp \sum_{r,r'}
 \frac{1}{(\varepsilon_{r'{\bf 