A rigorous, simple and efficient approach is derived in this paper for multipole expansion of a circularly symmetric Bessel beam. Different from the existing rigorous methods which are based on the plane wave spectrum of a Bessel beam, a straight-forward integral procedure is presented in a traditional way to obtain the analytical expressions of the expansion coefficients, also called beam shape coefficients (BSCs). The convergence and correctness of the BSCs are verified numerically in detail for both on-axis and off-axis cases. The results in this paper are useful in various analytical scattering theories, such as the generalized Lorenz-Mie theory and the Null-field method, when a Bessel beam is considered.