Abstract: The convolution quadrature based boundary element method (CQ-BEM) can preserve the kernel Dirac function when dealing with integral in time domain, forming a time convolution, which is a stable and effective numerical method in the time domain acoustic boundary element method. However, the traditional way to acquire the convolution integral coefficients in CQ-BEM requires large computation, long calculating time, and recalculation for different elements, which greatly reduces the efficiency of calculating the sound field in time domain. In order to solve the problem, this paper proposes the analytical expression and the numerical method of the coefficients based on polynomial theorem. The conversion theory of the coefficients between different elements is established so that the coefficients of different elements can be acquired in a single calculation, which greatly reduces the calculation amount. The sound field calculation efficiency of CQ-BEM method is then improved. In the numerical example of pulsating spherical source, the calculation time, relative error and growth rate of calculation time with the number of elements of traditional method and proposed method are compared. The results show that, under the same requirements, the calculation time of this method is more than 50% less than the traditional method, the relative error is more than 5 orders of magnitude less, and the growth rate of calculation time with the number of elements is only 2.34% of the traditional method, which can effectively improve the computational efficiency of CQ-BEM method in time domain sound field.