gives $$\begin{aligned}
 \sigma_{N}^{\rm indirect}(\omega) 
 &=g_s g_v
 \frac{e^2\hbar}{V} \sum_{\bf p} \sum_{r}
 \frac{\pi \delta (\hbar \omega - 2 (m_r +E_r))}{2(m_r +E_r)} 
 \frac{\varepsilon_{N+1-r{\bf p}}^{-}}{2E_r}
 \frac{\varepsilon_{r{\bf p}}^{+}}{2E_r} 
 (p_y)^2 \left( 
 \frac{1}{-m_r -E_r} - \frac{1}{m_r + E_r} \right)^2 \nn \\
 &=g_s g_v
 \frac{S}{V}
 \frac{e^2}{16\hbar} \sum_{r}\left(1- \frac{2m_r}{\hbar \omega -2m_r}\right)\Theta (\hbar \omega - 2 (m_r+|m_r|)).
 \label{eq:indirect}\end{aligned}$$ It looks as if $\sigma_{N}^{\rm indirect}(\omega)$ can be singular at $\hbar \omega = 2m_r$. However, due to the step function, $\sigma_{N}^{\rm indirect}(\omega)$ is a sum of (at most) discontinuous functions.

By combining Eqs. (\[eq:direct\]) and (\[eq:indirect\]), we obtain a simple formula of the bunched conductivity as $$\begin{aligned}
 \sigma_{N}(\omega) = \frac{S}{V} g_s g_v
 \frac{e^2}{16\hbar} \sum_{r=1}^N \left\{
 \left(1- \frac{2m_r}{\hbar \omega -2m_r}\right)\Theta (\hbar \omega - 2 (m_r+|m_r|))
 + \left( \frac{2m_r}{\hbar \omega} \right)^2 \Theta (\hbar \omega - 2 |m_r|) \right\}.
 \label{eq:bunchedc}\end{aligned}$$ Here, $V$ denotes the three dimensional volume of the system, $V=S$ for $N=1$ and $V=N S d$ for $N\ge 2$. This formula reproduces the results obtained previously for $N=2$ [@Abergel2007; @Koshino2008] and for small $N \le 10$. [@Min2009; @Orlita2010] The case of (AB stacking) graphite, i.e., $\lim_{N\to \infty} \sigma_N(\omega)$ can be analytically calculated and the result is shown in Appendix \[app:largeN\]. Our result differs slightly from the calculations obtained previously by several authors. [@Boyle1958; @Ichikawa1966]

