modes of electromagnetic fluctuations ($B_z=0$ and $E_z|_{\rho=R}=0$) which obey Dirichlet boundary condition are given by the Bessel functions (so that $E_z=J_\sigma(\alpha_{\sigma,n}\rho/R)\text{e}^{i\sigma\theta}\text{e}^{ikz}$), and the Bessel zeros specified by $\alpha_{\sigma,n}$, with $\sigma=0,\pm1,\pm2,...$ and $n=1,2,3,...$, naturally appear in the expression of the vacuum energy [@Jackson]. On the other hand, the TE modes ($E_z=0$ and $\frac{\partial B_z}{
\partial\rho}|_{\rho=R}=0$) inside the cylinder are compatible with the Neumann boundary condition, and they are given by $B_z=J_\sigma(\beta_{\sigma,n}\rho/R)\text{e}^{i\sigma\theta}\text{e}^{ikz}$ where $\beta_{\sigma,n}$s are positive zeros of the first derivatives of the Bessel functions of the first kind [@Jackson]. The two distinct class of modes, as shown in Fig. \[fig1\], contribute simultaneously to the vacuum energy of the conducting cylinder, which can be obtained from Eq. (\[eqn:11\]) to be [@Gosdzinsky; @Milton1999; @Jackson] $$\begin{aligned}
\label{eqn:12}
E= \frac{\hbar c L}{2}\sum_{\sigma,n} \int_{-\infty}^\infty \frac{dk}{2\pi}\bigg[\sqrt{\frac{\alpha_{\sigma,n}^2}{R^2}+k^2}+\sqrt{\frac{\beta_{\sigma,n}^2}{R^2}+k^2}\bigg],\end{aligned}$$ where the first and second terms in the square bracket represent contributions for the TM modes and TE modes respectively. Note that the integral over $k$ in the above expression has a logarithmic divergence which may not be easily removed by applying standard analytic continuation techniques. What work here are ultraviolet cutoff technique with Green’s dyadic formulation for the field strengths [@Milton1981], zeta regulation technique [@Gosdzinsky], mode-by-mode summation technique [@Milton1999], and, 